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: 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 : //# $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 0 : SpectralElement::SpectralElement(SpectralElement::Types type, const Vector<Double>& parms)
48 : : _type(type), _params(parms), _errors(parms.size(), 0),
49 0 : _fixed(parms.size(), false) {}
50 :
51 :
52 0 : SpectralElement::SpectralElement(const SpectralElement &other)
53 0 : : _type(other._type), _params(other._params.copy()), _errors(other._errors.copy()),
54 0 : _fixed(other._fixed.copy()),
55 : _function(
56 : std::shared_ptr<Function<Double, Double> >(
57 0 : other._function->clone()
58 : )
59 0 : ) {}
60 :
61 0 : SpectralElement::~SpectralElement() {}
62 :
63 0 : SpectralElement &SpectralElement::operator=(
64 : const SpectralElement &other
65 : ) {
66 0 : if (this != &other) {
67 0 : _type = other._type;
68 0 : uInt n = other._params.size();
69 0 : _params.resize(n);
70 0 : _params = other._params.copy();
71 0 : _errors.resize(n);
72 0 : _errors = other._errors.copy();
73 0 : _fixed.resize(n);
74 0 : _fixed = other._fixed.copy();
75 0 : _function = std::shared_ptr<Function<Double, Double> >(
76 0 : other._function->clone()
77 0 : );
78 : }
79 0 : 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 0 : 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 0 : };
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 0 : nall = SpectralElement::N_Types;
128 0 : typ = oname;
129 0 : return tname;
130 : }
131 :
132 0 : void SpectralElement::_set(const Vector<Double>& params) {
133 0 : _params = params.copy();
134 0 : for (uInt i=0; i<params.size(); i++) {
135 0 : (*_function)[i] = params[i];
136 : }
137 0 : }
138 :
139 0 : void SpectralElement::_setType(const SpectralElement::Types type) {
140 0 : _type = type;
141 0 : }
142 :
143 0 : const String &SpectralElement::fromType(SpectralElement::Types tp) {
144 : Int nall;
145 : const SpectralElement::Types *typ;
146 0 : const String *const tname = SpectralElement::allTypes(nall, typ);
147 0 : 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 0 : void SpectralElement::_setFunction(
168 : const std::shared_ptr<Function<Double, Double> >& f
169 : ) {
170 0 : _function = f;
171 0 : }
172 :
173 0 : Double SpectralElement::operator()(const Double x) const {
174 0 : return (*_function)(x);
175 : }
176 :
177 0 : void SpectralElement::get(Vector<Double> ¶m) const {
178 0 : param = _params.copy();
179 0 : }
180 :
181 0 : Vector<Double> SpectralElement::get() const {
182 0 : return _params.copy();
183 : }
184 :
185 0 : void SpectralElement::getError(Vector<Double> &err) const {
186 0 : err = _errors.copy();
187 0 : }
188 :
189 0 : Vector<Double> SpectralElement::getError() const {
190 0 : return _errors.copy();
191 : }
192 :
193 0 : void SpectralElement::setError(const Vector<Double> &err) {
194 0 : if (err.nelements() != _params.nelements()) {
195 : throw(
196 : AipsError(
197 : "SpectralElement: setting incorrect "
198 : "number of errors in the element"
199 : )
200 0 : );
201 : }
202 0 : _errors = err.copy();
203 0 : }
204 :
205 0 : void SpectralElement::fix(const Vector<Bool> &fix) {
206 0 : if (fix.nelements() != _params.nelements()) {
207 : throw(
208 : AipsError(
209 : "SpectralElement: setting incorrect number of fixed "
210 : "in the element"
211 : )
212 0 : );
213 : }
214 0 : _fixed = fix.copy();
215 0 : for (uInt i=0; i<_fixed.size(); i++) {
216 0 : _function->mask(i) = fix[i];
217 : }
218 0 : }
219 :
220 0 : const Vector<Bool>& SpectralElement::fixed() const {
221 0 : 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 : }
233 :
234 0 : void SpectralElement::set(const Vector<Double>& params) {
235 0 : if (params.size() != get().size()) {
236 : throw AipsError(
237 : "Input number of parameters does not match "
238 : "the current number of parameters"
239 0 : );
240 : }
241 0 : _set(params);
242 0 : }
243 :
244 : } //# NAMESPACE CASA - END
245 :
246 :
|
