Line data    Source code

       1             : //# Spectral2Element.cc: Record conversion for SpectralElement
       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: Spectral2Element.cc 20652 2009-07-06 05:04:32Z Malte.Marquarding $
      27             : 
      28             : //# Includes
      29             : 
      30             : #include <casacore/casa/Containers/Record.h>
      31             : #include <memory>
      32             : #include <components/SpectralComponents/CompiledSpectralElement.h>
      33             : #include <components/SpectralComponents/GaussianSpectralElement.h>
      34             : #include <components/SpectralComponents/GaussianMultipletSpectralElement.h>
      35             : #include <components/SpectralComponents/LogTransformedPolynomialSpectralElement.h>
      36             : #include <components/SpectralComponents/LorentzianSpectralElement.h>
      37             : #include <components/SpectralComponents/PolynomialSpectralElement.h>
      38             : #include <components/SpectralComponents/PowerLogPolynomialSpectralElement.h>
      39             : #include <components/SpectralComponents/SpectralElementFactory.h>
      40             : 
      41             : using namespace casacore;
      42             : namespace casa { //# NAMESPACE CASA - BEGIN
      43             : 
      44           0 : SpectralElement* SpectralElementFactory::fromRecord(
      45             :         const RecordInterface &in
      46             : ) {
      47           0 :         std::unique_ptr<SpectralElement> specEl;
      48           0 :         String origin = "SpectralElementFactory::fromRecord: ";
      49           0 :         if (
      50           0 :                 ! in.isDefined("type")
      51           0 :                 || in.type(in.idToNumber(RecordFieldId("type"))) != TpString
      52             :         ) {
      53           0 :                 throw AipsError("Record does not represent a SpectralElement");
      54             :         }
      55           0 :         String stp;
      56             :         SpectralElement::Types tp;
      57           0 :         in.get(RecordFieldId("type"), stp);
      58           0 :         if (!SpectralElement::toType(tp, stp)) {
      59           0 :                 throw AipsError("Unknown spectral type in SpectralElement::fromRecord\n");
      60             :         }
      61             : 
      62           0 :         Vector<Double> errs;
      63             : 
      64             :         // Get the errors if defined in record
      65             : 
      66           0 :         if (in.isDefined("errors")) {
      67           0 :                 DataType dataType = in.dataType("errors");
      68           0 :                 if (dataType == TpArrayDouble) {
      69           0 :                         in.get("errors", errs);
      70             :                 }
      71           0 :                 else if (dataType == TpArrayFloat) {
      72           0 :                         Vector<Float> v;
      73           0 :                         in.get("errors", v);
      74           0 :                         errs.resize(v.nelements());
      75           0 :                         convertArray(errs, v);
      76             :                 }
      77           0 :                 else if (dataType == TpArrayInt) {
      78           0 :                         Vector<Int> v;
      79           0 :                         in.get("errors", v);
      80           0 :                         errs.resize(v.nelements());
      81           0 :                         convertArray(errs, v);
      82             :                 }
      83             :                 else {
      84             :                         throw AipsError(
      85             :                                 "SpectralElement::fromRecord: errors field "
      86             :                                 "must be double, float or int\n"
      87           0 :                         );
      88             :                 }
      89             :         }
      90             : 
      91           0 :         Vector<Double> param;
      92           0 :         if (in.isDefined(("parameters"))) {
      93           0 :                 DataType dataType = in.dataType("parameters");
      94           0 :                 if (dataType == TpArrayDouble) {
      95           0 :                         in.get("parameters", param);
      96             :                 }
      97           0 :                 else if (dataType == TpArrayFloat) {
      98           0 :                         Vector<Float> v;
      99           0 :                         in.get("parameters", v);
     100           0 :                         param.resize(v.nelements());
     101           0 :                         convertArray(param, v);
     102             :                 }
     103           0 :                 else if (dataType == TpArrayInt) {
     104           0 :                         Vector<Int> v;
     105           0 :                         in.get("parameters", v);
     106           0 :                         param.resize(v.nelements());
     107           0 :                         convertArray(param, v);
     108             :                 }
     109             :                 else {
     110             :                         throw AipsError(
     111           0 :                                 origin +
     112             :                                 "SpectralElement::fromRecord: parameters field "
     113             :                                 "must be double, float or int\n"
     114           0 :                         );
     115             :                 }
     116             :         }
     117             : 
     118             :         // Make sizes of errors and params equal
     119           0 :         if (errs.nelements() == 0) {
     120           0 :                 errs.resize(param.nelements());
     121           0 :                 errs = 0.0;
     122             :         }
     123           0 :         if (errs.nelements() != param.nelements()) {
     124             :                 throw AipsError(
     125           0 :                         origin +
     126             :                         "SpectralElement::fromRecord must have equal lengths "
     127             :                         "for parameters and errors fields"
     128           0 :                 );
     129             :         }
     130           0 :         switch (tp) {
     131           0 :         case SpectralElement::GAUSSIAN:
     132           0 :                 if (param.nelements() != 3) {
     133             :                         throw AipsError(
     134           0 :                                 origin + "Illegal number of parameters for Gaussian element"
     135           0 :                         );
     136             :                 }
     137           0 :                 if (param(2) <= 0.0) {
     138             :                         throw AipsError(
     139           0 :                                 origin +
     140             :                                 "The width of a Gaussian element must be positive"
     141           0 :                         );
     142             :                 }
     143           0 :                 param(2) = GaussianSpectralElement::sigmaFromFWHM (param(2));
     144           0 :                 errs(2) = GaussianSpectralElement::sigmaFromFWHM (errs(2));
     145           0 :                 specEl = std::unique_ptr<SpectralElement>(new GaussianSpectralElement(param(0), param(1), param(2)));
     146           0 :                 specEl->setError(errs);
     147           0 :                 break;
     148           0 :         case SpectralElement::LORENTZIAN:
     149           0 :                 if (param.nelements() != 3) {
     150             :                         throw AipsError(
     151             :                                 "Illegal number of parameters for Lorentzian element"
     152           0 :                         );
     153             :                 }
     154           0 :                 if (param(2) <= 0.0) {
     155             :                         throw AipsError(
     156             :                                 "The width of a Lorentzian element must be positive"
     157           0 :                         );
     158             :                 }
     159           0 :                 specEl = std::unique_ptr<SpectralElement>(new LorentzianSpectralElement(param(0), param(1), param(2)));
     160           0 :                 specEl->setError(errs);
     161           0 :                 break;
     162           0 :         case SpectralElement::POLYNOMIAL:
     163           0 :                 if (param.nelements() == 0) {
     164             :                         throw AipsError(
     165             :                                 "Polynomial spectral element must have order "
     166             :                                 "of at least zero"
     167           0 :                         );
     168             :                 }
     169           0 :                 specEl = std::unique_ptr<SpectralElement>(new PolynomialSpectralElement(param.nelements() - 1));
     170           0 :                 specEl->set(param);
     171           0 :                 specEl->setError(errs);
     172           0 :                 break;
     173           0 :         case SpectralElement::COMPILED:
     174           0 :                 if (
     175           0 :                                 in.isDefined("compiled")
     176           0 :                                 && in.type(in.idToNumber(RecordFieldId("compiled"))) == TpString
     177             :                 ) {
     178           0 :                         String function;
     179           0 :                         in.get(RecordFieldId("compiled"), function);
     180           0 :                         specEl = std::unique_ptr<SpectralElement>(new CompiledSpectralElement(function, param));
     181           0 :                         specEl->setError(errs);
     182             :                 }
     183             :                 else {
     184             :                         throw AipsError(
     185             :                                 "No compiled string in SpectralElement::fromRecord\n"
     186           0 :                         );
     187             :                 }
     188           0 :                 break;
     189             : 
     190           0 :         case SpectralElement::GMULTIPLET:
     191             :         {
     192           0 :                 if (! in.isDefined("gaussians")) {
     193           0 :                         throw AipsError("gaussians not defined in record");
     194             :                 }
     195           0 :                 if (! in.isDefined("fixedMatrix")) {
     196           0 :                         throw AipsError("fixed matrix not defined in record");
     197             :                 }
     198           0 :                 Record gaussians = in.asRecord("gaussians");
     199           0 :                 uInt i = 0;
     200           0 :                 std::vector<GaussianSpectralElement> comps(0);
     201             :                 while(true) {
     202           0 :                         String id = "*" + String::toString(i);
     203           0 :                         if (gaussians.isDefined(id)) {
     204           0 :                                 std::unique_ptr<SpectralElement> gauss(fromRecord(gaussians.asRecord(id)));
     205           0 :                                 comps.push_back(
     206             :                                         *dynamic_cast<GaussianSpectralElement*>(
     207           0 :                                                 gauss.get()
     208           0 :                                         )
     209             :                                 );
     210           0 :                                 i++;
     211             :                         }
     212             :                         else {
     213           0 :                                 break;
     214             :                         }
     215           0 :                 }
     216           0 :                 Matrix<Double> fixedMatrix = in.asArrayDouble("fixedMatrix");
     217           0 :                 fixedMatrix.reform(IPosition(2, comps.size()-1, 3));
     218           0 :                 specEl = std::unique_ptr<SpectralElement>(new GaussianMultipletSpectralElement(comps, fixedMatrix));
     219             :         }
     220           0 :         break;
     221             : 
     222           0 :     case SpectralElement::POWERLOGPOLY: {
     223           0 :                 specEl = std::unique_ptr<SpectralElement>(new PowerLogPolynomialSpectralElement(param));
     224           0 :                 specEl->set(param);
     225           0 :                 specEl->setError(errs);
     226             :         }
     227           0 :         break;
     228             : 
     229           0 :     case SpectralElement::LOGTRANSPOLY: {
     230           0 :             specEl = std::unique_ptr<SpectralElement>(new LogTransformedPolynomialSpectralElement(param));
     231           0 :                 specEl->setError(errs);
     232             :         }
     233           0 :         break;
     234             : 
     235           0 :         default:
     236             :                 throw AipsError(
     237             :                         "Unhandled or illegal spectral element record in "
     238             :                         "SpectralElementFactory::fromRecord\n"
     239           0 :                 );
     240             :         }
     241             : 
     242           0 :         if (in.isDefined("fixed")) {
     243           0 :                 specEl->fix(in.asArrayBool("fixed"));
     244             :         }
     245             :     // ready to return, fish out the pointer and return it without deleting it
     246           0 :     SpectralElement *sp = specEl.get();
     247           0 :     specEl.release( );
     248           0 :         return sp;
     249             : }
     250             : 
     251             : 
     252             : } //# NAMESPACE CASA - END
     253             :