Line data Source code
1 : //# Copyright (C) 1993,1994,1995,1996,1999,2001
2 : //# Associated Universities, Inc. Washington DC, USA.
3 : //#
4 : //# This library is free software; you can redistribute it and/or modify it
5 : //# under the terms of the GNU Library General Public License as published by
6 : //# the Free Software Foundation; either version 2 of the License, or (at your
7 : //# option) any later version.
8 : //#
9 : //# This library is distributed in the hope that it will be useful, but WITHOUT
10 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
12 : //# License for more details.
13 : //#
14 : //# You should have received a copy of the GNU Library General Public License
15 : //# along with this library; if not, write to the Free Software Foundation,
16 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
17 : //#
18 : //# Correspondence concerning AIPS++ should be addressed as follows:
19 : //# Internet email: aips2-request@nrao.edu.
20 : //# Postal address: AIPS++ Project Office
21 : //# National Radio Astronomy Observatory
22 : //# 520 Edgemont Road
23 : //# Charlottesville, VA 22903-2475 USA
24 : //#
25 :
26 : #include <components/SpectralComponents/SpectralListFactory.h>
27 :
28 : #include <components/SpectralComponents/GaussianSpectralElement.h>
29 : #include <components/SpectralComponents/GaussianMultipletSpectralElement.h>
30 : #include <components/SpectralComponents/LorentzianSpectralElement.h>
31 : #include <components/SpectralComponents/PowerLogPolynomialSpectralElement.h>
32 : #include <components/SpectralComponents/LogTransformedPolynomialSpectralElement.h>
33 :
34 : #include <stdcasa/StdCasa/CasacSupport.h>
35 :
36 : using namespace casacore;
37 : using namespace casac;
38 :
39 : using namespace casacore;
40 : namespace casa {
41 :
42 0 : SpectralList SpectralListFactory::create(
43 : LogIO& log, const variant& pampest,
44 : const variant& pcenterest, const variant& pfwhmest,
45 : const variant& pfix, const variant& gmncomps,
46 : const variant& gmampcon, const variant& gmcentercon,
47 : const variant& gmfwhmcon, const vector<double>& gmampest,
48 : const vector<double>& gmcenterest, const vector<double>& gmfwhmest,
49 : const variant& gmfix, const variant& pfunc,
50 : const variant& plpest, const variant& plpfix,
51 : const variant& ltpest, const variant& ltpfix
52 : ) {
53 0 : vector<double> myampest = toVectorDouble(pampest, "pampest");
54 0 : vector<double> mycenterest = toVectorDouble(pcenterest, "pcenterest");
55 0 : vector<double> myfwhmest = toVectorDouble(pfwhmest, "pfwhmest");
56 0 : vector<string> myfix = toVectorString(pfix, "pfix");
57 0 : vector<string> myfunc = toVectorString(pfunc, "pfunc");
58 0 : vector<double> myplpest = toVectorDouble(plpest, "plpest");
59 0 : vector<bool> myplpfix = toVectorBool(plpfix, "plpfix");
60 0 : vector<double> myltpest = toVectorDouble(ltpest, "ltpest");
61 0 : vector<bool> myltpfix = toVectorBool(ltpfix, "ltpfix");
62 :
63 0 : vector<int> mygmncomps = gmncomps.type() == variant::INT
64 0 : && gmncomps.toInt() == 0
65 0 : ? vector<int>(0)
66 0 : : toVectorInt(gmncomps, "gmncomps");
67 0 : vector<double> mygmampcon = toVectorDouble(gmampcon, "gmampcon");
68 0 : vector<double> mygmcentercon = toVectorDouble(gmcentercon, "gmcentercon");
69 0 : vector<double> mygmfwhmcon = toVectorDouble(gmfwhmcon, "gmfwhmcon");
70 0 : vector<string> mygmfix = toVectorString(gmfix, "gmfix");
71 : Bool makeSpectralList = (
72 0 : mygmncomps.size() + myplpest.size() + myltpest.size() > 0
73 0 : || ! (
74 0 : myampest.size() == 0
75 0 : && mycenterest.size() == 0
76 0 : && myfwhmest.size() == 0
77 : )
78 0 : );
79 0 : SpectralList spectralList;
80 0 : Bool gfixSpecified = (myfix.size() > 0);
81 0 : if (! makeSpectralList) {
82 0 : if (gfixSpecified) {
83 : log << "The fix array is specified but no corresponding estimates are "
84 : << "set via ampest, centerest, and fwhmest. The fix array will be ignored."
85 0 : << LogIO::WARN;
86 : }
87 0 : return spectralList;
88 : }
89 0 : uInt mynpcf = myampest.size();
90 0 : if (myfunc.size() == 0) {
91 0 : myfunc.resize(myampest.size(), "G");
92 : }
93 0 : ThrowIf(
94 : mycenterest.size() != mynpcf
95 : || myfwhmest.size() != mynpcf
96 : || myfunc.size() != mynpcf,
97 : "pampest, pcenterest, pfwhmest, and pfunc arrays "
98 : "must all be the same length"
99 : );
100 0 : ThrowIf(
101 : gfixSpecified && myfix.size() != mynpcf,
102 : "If the gfix array is specified the number of "
103 : "elements it has must be the same as the number of elements "
104 : "in the ampest array even if some elements are empty strings"
105 : );
106 0 : if (mygmncomps.size() > 0) {
107 0 : _addGaussianMultiplets(
108 : spectralList, log, mygmncomps, mygmampcon,
109 : mygmcentercon, mygmfwhmcon,
110 : gmampest, gmcenterest, gmfwhmest,
111 : mygmfix
112 : );
113 : }
114 0 : for (uInt i=0; i<mynpcf; i++) {
115 0 : String func(myfunc[i]);
116 0 : func.upcase();
117 0 : Bool doGauss = func.startsWith("G");
118 0 : Bool doLorentz = func.startsWith("L");
119 0 : if (! doGauss && ! doLorentz) {
120 0 : log << myfunc[i] << " does not minimally match 'gaussian' or 'lorentzian'"
121 0 : << LogIO::EXCEPTION;
122 : }
123 : std::unique_ptr<PCFSpectralElement> pcf(
124 : doGauss
125 : ? dynamic_cast<PCFSpectralElement*>(
126 : new GaussianSpectralElement(
127 0 : myampest[i], mycenterest[i],
128 0 : GaussianSpectralElement::sigmaFromFWHM(myfwhmest[i])
129 0 : )
130 : )
131 : : doLorentz
132 0 : ? dynamic_cast<PCFSpectralElement*>(
133 : new LorentzianSpectralElement(
134 0 : myampest[i], mycenterest[i], myfwhmest[i]
135 0 : )
136 : )
137 : : 0
138 0 : );
139 0 : if (gfixSpecified) {
140 0 : pcf->fixByString(myfix[i]);
141 : }
142 0 : if (! spectralList.add(*pcf)) {
143 : log << "Unable to add element to spectral list"
144 0 : << LogIO::EXCEPTION;
145 : }
146 : }
147 0 : if (myplpest.size() > 0) {
148 0 : _addPowerLogPolynomial(
149 : spectralList, log, myplpest, myplpfix
150 : );
151 : }
152 0 : if (myltpest.size() > 0) {
153 0 : _addLogTransformedPolynomial(
154 : spectralList, myltpest, myltpfix
155 : );
156 : }
157 0 : return spectralList;
158 : }
159 :
160 0 : void SpectralListFactory::_addGaussianMultiplets(
161 : SpectralList& spectralList,
162 : LogIO& log,
163 : const vector<int>& mygmncomps,
164 : vector<double>& mygmampcon,
165 : vector<double>& mygmcentercon,
166 : vector<double>& mygmfwhmcon,
167 : const vector<double>& gmampest,
168 : const vector<double>& gmcenterest, const vector<double>& gmfwhmest,
169 : const vector<string>& mygmfix
170 : ) {
171 0 : uInt sum = 0;
172 0 : for (uInt i=0; i<mygmncomps.size(); i++) {
173 0 : if (mygmncomps[i] <= 0) {
174 0 : log << "All elements of gmncomps must be greater than 0" << LogIO::EXCEPTION;
175 : }
176 0 : sum += mygmncomps[i];
177 : }
178 0 : if (gmampest.size() != sum) {
179 : log << "gmampest must have exactly "
180 0 : << sum << " elements" << LogIO::EXCEPTION;
181 : }
182 0 : if (gmcenterest.size() != sum) {
183 : log << "gmcenterest must have exactly "
184 0 : << sum << " elements" << LogIO::EXCEPTION;
185 : }
186 0 : if (gmfwhmest.size() != sum) {
187 : log << "gmfwhmest must have exactly "
188 0 : << sum << " elements" << LogIO::EXCEPTION;
189 : }
190 0 : if (mygmfix.size() > 0 && mygmfix.size() != sum) {
191 : log << "gmfwhmest must have either zero or " << sum
192 : << " elements, even if some are empty strings."
193 0 : << LogIO::EXCEPTION;
194 : }
195 0 : uInt nConstraints = sum - mygmncomps.size();
196 0 : if (mygmampcon.size() == 0) {
197 0 : mygmampcon = vector<double>(nConstraints, 0);
198 : }
199 0 : else if (
200 0 : mygmampcon.size() > 0
201 0 : && mygmampcon.size() != nConstraints
202 : ) {
203 : log << "If specified, gmampcon must have exactly "
204 : << nConstraints << " elements, even if some are zero"
205 0 : << LogIO::EXCEPTION;
206 : }
207 0 : if (mygmcentercon.size() == 0) {
208 0 : mygmcentercon = vector<double>(nConstraints, 0);
209 : }
210 0 : else if (
211 0 : mygmcentercon.size() > 0
212 0 : && mygmcentercon.size() != nConstraints
213 : ) {
214 : log << "If specified, gmcentercon must have exactly "
215 : << nConstraints << " elements, even if some are zero"
216 0 : << LogIO::EXCEPTION;
217 : }
218 0 : if (mygmfwhmcon.size() == 0) {
219 0 : mygmfwhmcon = vector<double>(nConstraints, 0);
220 : }
221 0 : else if (
222 0 : mygmfwhmcon.size() > 0
223 0 : && mygmfwhmcon.size() != nConstraints
224 : ) {
225 : log << "If specified, gmfwhmcon must have exactly "
226 : << nConstraints << " elements, even if some are zero"
227 0 : << LogIO::EXCEPTION;
228 : }
229 0 : uInt compNumber = 0;
230 0 : for (uInt i=0; i<mygmncomps.size(); i++) {
231 0 : if (mygmncomps[i] < 0) {
232 : log << "All elements of gmncomps must be positive"
233 0 : << LogIO::EXCEPTION;
234 : }
235 0 : vector<GaussianSpectralElement> g(mygmncomps[i]);
236 0 : Matrix<Double> constraints(mygmncomps[i] - 1, 3);
237 0 : for (uInt j=0; j<(uInt)mygmncomps[i]; j++) {
238 0 : g[j] = GaussianSpectralElement(
239 0 : gmampest[compNumber], gmcenterest[compNumber],
240 0 : GaussianSpectralElement::sigmaFromFWHM(gmfwhmest[compNumber])
241 0 : );
242 0 : if (mygmfix.size() > 0) {
243 0 : g[j].fixByString(mygmfix[compNumber]);
244 : }
245 0 : if (j > 0) {
246 0 : constraints(j-1, 0) = mygmampcon[compNumber - (i+1)];
247 0 : constraints(j-1, 1) = mygmcentercon[compNumber - (i+1)];
248 0 : constraints(j-1, 2) = mygmfwhmcon[compNumber - (i+1)];
249 : }
250 0 : compNumber++;
251 : }
252 0 : if (
253 0 : ! spectralList.add(
254 0 : GaussianMultipletSpectralElement(g, constraints)
255 : )
256 : ) {
257 : log << "Unable to add gaussian multiplet to spectral list"
258 0 : << LogIO::EXCEPTION;
259 : }
260 : }
261 0 : }
262 :
263 0 : void SpectralListFactory::_addPowerLogPolynomial(
264 : SpectralList& spectralList,
265 : LogIO& log, vector<double>& myplpest,
266 : vector<bool>& myplpfix
267 : ) {
268 0 : uInt nest = myplpest.size();
269 0 : if (nest < 2) {
270 : log << "Number of elements in the power log polynomial estimates list must be at least 2"
271 0 : << LogIO::EXCEPTION;
272 : }
273 0 : uInt nfix = myplpfix.size();
274 0 : if (nfix == 0) {
275 0 : myplpfix = vector<bool>(nest, false);
276 : }
277 0 : else if (nfix != myplpest.size()) {
278 : log << "Number of elements in the power logarithmic polynomial fixed parameter list must "
279 : << "either be 0 or equal to the number of elements in the estimates list"
280 0 : << LogIO::EXCEPTION;
281 : }
282 0 : Vector<double> const myplpestV(myplpest);
283 0 : Vector<bool> const myplpfixV(myplpfix);
284 0 : PowerLogPolynomialSpectralElement plp(myplpestV);
285 0 : plp.fix(myplpfixV);
286 0 : spectralList.add(plp);
287 0 : }
288 :
289 0 : void SpectralListFactory::_addLogTransformedPolynomial(
290 : SpectralList& spectralList,
291 : vector<double>& myltpest,
292 : vector<bool>& myltpfix
293 : ) {
294 0 : uInt nest = myltpest.size();
295 0 : ThrowIf(
296 : nest < 2,
297 : "Number of elements in the log transformed polynomial "
298 : "estimates list must be at least 2"
299 : );
300 0 : uInt nfix = myltpfix.size();
301 0 : if (nfix == 0) {
302 0 : myltpfix = vector<bool>(nest, false);
303 : }
304 : else {
305 0 : ThrowIf(
306 : nfix != myltpest.size(),
307 : "Number of elements in the logarithmic transformed polynomial "
308 : "fixed parameter list must either be 0 or equal to the number "
309 : "of elements in the estimates list"
310 : );
311 : }
312 0 : Vector<double> const myltpestV(myltpest);
313 0 : Vector<bool> const myltpfixV(myltpfix);
314 0 : LogTransformedPolynomialSpectralElement ltp(myltpestV);
315 0 : ltp.fix(myltpfixV);
316 0 : spectralList.add(ltp);
317 0 : }
318 :
319 : using namespace casacore;
320 : } // end namespace casa
|
