jacobian.h

Go to the documentation of this file.

/* Copyright (C) 2004-2012 Mattias Ekstrom <ekstrom@rss.chalmers.se>
 
   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.
 
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
 
   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */
 
#ifndef jacobian_h
#define jacobian_h
 
#include "species_tags.h"
#include "agenda_class.h"
#include "array.h"
#include "bifstream.h"
#include "enums.h"
#include "interpolation.h"
#include "logic.h"
#include "matpackI.h"
#include "methods.h"
#include "mystring.h"
#include "ppath.h"
#include "quantum.h"  
#include <iostream>
#include <map>
#include <stdexcept>
 
namespace Jacobian {
 
ENUMCLASS(Type, char,
          Atm,
          Line,
          Sensor,
          Special
          )
 
 
ENUMCLASS(Atm, char,
          Temperature,
          WindMagnitude, WindU, WindV, WindW,
          MagneticMagnitude, MagneticU, MagneticV, MagneticW,
          Electrons,
          Particulates
          )
 
ENUMCLASS(Line, char,
          VMR,
          Strength,
          Center,
          ShapeG0X0, ShapeG0X1, ShapeG0X2, ShapeG0X3,
          ShapeD0X0, ShapeD0X1, ShapeD0X2, ShapeD0X3,
          ShapeG2X0, ShapeG2X1, ShapeG2X2, ShapeG2X3,
          ShapeD2X0, ShapeD2X1, ShapeD2X2, ShapeD2X3,
          ShapeFVCX0, ShapeFVCX1, ShapeFVCX2, ShapeFVCX3,
          ShapeETAX0, ShapeETAX1, ShapeETAX2, ShapeETAX3,
          ShapeYX0, ShapeYX1, ShapeYX2, ShapeYX3,
          ShapeGX0, ShapeGX1, ShapeGX2, ShapeGX3,
          ShapeDVX0, ShapeDVX1, ShapeDVX2, ShapeDVX3,
          NLTE,
          ECS_A, ECS_B, ECS_GAMMA, ECS_DC
          )
 
ENUMCLASS(Sensor, char,
          FrequencyShift,
          FrequencyStretch,
          Polyfit,
          Sinefit,
          PointingZenithInterp,
          PointingZenithRecalc
          )
 
ENUMCLASS(Special, char,
          ArrayOfSpeciesTagVMR,
          ScatteringString,
          SurfaceString
          )
 
class Target {
private:
  Type mtype{Type::FINAL};
  
  Atm matm{Atm::FINAL};
  
  Line mline{Line::FINAL};
  
  Sensor msensor{Sensor::FINAL};
  
  Special mspecial{Special::FINAL};
  
  Numeric mperturbation{std::numeric_limits<Numeric>::quiet_NaN()};
  
  QuantumIdentifier mqid{};
  
  ArrayOfSpeciesTag maostid{0};
  
  String msid{};
  
  Species::Species mspecid{Species::Species::FINAL};
public:
  explicit Target (Atm type) :
  mtype(Type::Atm), matm(type) {
    ARTS_ASSERT(good_enum(type))
  }
  
  explicit Target (Line type, const QuantumIdentifier& qid, Species::Species specid) :
  mtype(Type::Line), mline(type), mqid(qid), mspecid(specid) {
    ARTS_ASSERT(good_enum(type))
  }
  
  explicit Target (Sensor type) :
  mtype(Type::Sensor), msensor(type) {
    ARTS_ASSERT(good_enum(type))
  }
  
  explicit Target (Special type, const ArrayOfSpeciesTag& aostid) :
  mtype(Type::Special), mspecial(type), maostid(aostid) {
    ARTS_ASSERT (type == Special::ArrayOfSpeciesTagVMR,
      "Only for Special::ArrayOfSpeciesTagVMR, but you fed: ", mspecial)
  }
  
  explicit Target (Special type, const String& sid) :
  mtype(Type::Special), mspecial(type), msid(sid) {
    ARTS_ASSERT (type == Special::SurfaceString or type == Special::ScatteringString,
      "Only for Special::ScatteringString or Special::SurfaceString, but you fed: ", mspecial)
  }
  
  explicit Target () = default;
  
  void Perturbation(Numeric x) noexcept {mperturbation=x;}
  Numeric& Perturbation() noexcept {return mperturbation;}
  [[nodiscard]] Numeric Perturbation() const noexcept {return mperturbation;}
  
  void QuantumIdentity(const QuantumIdentifier& x) noexcept {mqid=x;}
  QuantumIdentifier& QuantumIdentity() noexcept {return mqid;}
  [[nodiscard]] const QuantumIdentifier& QuantumIdentity() const noexcept {return mqid;}
  
  void SpeciesList(const ArrayOfSpeciesTag& x) noexcept {maostid = x;}
  ArrayOfSpeciesTag& SpeciesList() noexcept {return maostid;}
  [[nodiscard]] const ArrayOfSpeciesTag& SpeciesList() const noexcept {return maostid;}
  
  void StringKey(const String& x) noexcept {msid = x;}
  String& StringKey() noexcept {return msid;}
  const String& StringKey() const noexcept {return msid;}
  
  [[nodiscard]] Atm AtmType() const noexcept {return matm;}
  
  [[nodiscard]] Line LineType() const noexcept {return mline;}
  
  [[nodiscard]] Sensor SensorType() const noexcept {return msensor;}
  
  [[nodiscard]] Special SpecialType() const noexcept {return mspecial;}
  
  bool operator==(Atm other) const noexcept {return other == matm;}
  
  bool operator==(Line other) const noexcept {return other == mline;}
  
  bool operator==(Sensor other) const noexcept {return other == msensor;}
  
  bool operator==(Special other) const noexcept {return other == mspecial;}
  
  bool operator==(Type other) const noexcept {return other == mtype;}
  
  [[nodiscard]] bool sameTargetType(const Target& other) const noexcept {
    return mtype    == other.mtype    and
           matm     == other.matm     and
           mline    == other.mline    and
           msensor  == other.msensor  and
           mspecial == other.mspecial;
  }
  
  std::string_view TargetType() const noexcept {return toString(mtype);}
  
  void TargetType(const std::string_view& s) noexcept {mtype = toType(s);}
  
  void TargetSubType(const std::string_view& s) noexcept {
    matm = Atm::FINAL;
    mline = Line::FINAL;
    msensor = Sensor::FINAL;
    mspecial = Special::FINAL;
    
    switch (mtype) {
      case Type::Atm: matm = toAtm(s); break;
      case Type::Line: mline = toLine(s); break;
      case Type::Sensor: msensor = toSensor(s); break;
      case Type::Special: mspecial = toSpecial(s); break;
      case Type::FINAL: { /* leave last, don't use default */ }
    }
  }
  
  std::string_view TargetSubType() const noexcept {
    switch (mtype) {
      case Type::Atm: return toString(matm);
      case Type::Line: return toString(mline);
      case Type::Sensor: return toString(msensor);
      case Type::Special: return toString(mspecial);
      case Type::FINAL: { /* leave last, don't use default */ }
    }
    return "BAD SUBTYPE";
  }
  
  bool TargetTypeOK() const noexcept {return good_enum(mtype);}
 
  bool TargetSubTypeOK() const noexcept {
    // We can only hold one valid enum at a time, and it must be of the correct type
    if (1 == (good_enum(mspecial) + good_enum(msensor) + good_enum(mline) + good_enum(matm))) {
      switch (mtype) {
        case Type::Special:
          return good_enum(mspecial);
        case Type::Sensor:
          return good_enum(msensor);
        case Type::Line:
          return good_enum(mline);
        case Type::Atm:
          return good_enum(matm);
        case Type::FINAL: { /* leave last, don't use default */ }
      }
    }
    
    return false;
  }
  
  bool isSpeciesVMR() const noexcept {
    return mline == Line::VMR or mspecial == Special::ArrayOfSpeciesTagVMR;
  }
  
  bool isWind() const noexcept {
    return matm == Atm::WindMagnitude or 
           matm == Atm::WindU or
           matm == Atm::WindV or
           matm == Atm::WindW;
  }
  
  bool isMagnetic() const noexcept {
    return matm == Atm::MagneticMagnitude or 
           matm == Atm::MagneticU or
           matm == Atm::MagneticV or
           matm == Atm::MagneticW;
  }
  
  bool isFrequency() const noexcept {
    return msensor == Sensor::FrequencyStretch or 
           msensor == Sensor::FrequencyShift;
  }
  
  bool isPointing() const noexcept {
    return msensor == Sensor::PointingZenithInterp or 
           msensor == Sensor::PointingZenithRecalc;
  }
  
  bool needQuantumIdentity() const noexcept {
    return mtype == Type::Line;
  }
  
  bool needArrayOfSpeciesTag() const noexcept {
    return mspecial == Special::ArrayOfSpeciesTagVMR;
  }
  
  bool needString() const noexcept {
    return mspecial == Special::ScatteringString or
           mspecial == Special::SurfaceString;
  }
  
  const Species::Species& LineSpecies() const noexcept {return mspecid;}
  Species::Species& LineSpecies() noexcept {return mspecid;}
  void LineSpecies(Species::Species x) noexcept {mspecid = x;}
};  // Target
 
std::ostream& operator<<(std::ostream& os, const Target& x);
};  // Jacobian
using JacobianTarget = Jacobian::Target;
using ArrayOfJacobianTarget = Array<Jacobian::Target>;
 
 
class RetrievalQuantity {
 public:
  RetrievalQuantity()
      : msubtag(),
        msubsubtag(),
        mmode(),
        mgrids(),
        mjac() { /* Nothing to do here. */
  }
 
  RetrievalQuantity(const Jacobian::Target& target,
                    const String& subtag,
                    const String& subsubtag,
                    const String& mode,
                    const Numeric& perturbation,
                    const ArrayOfVector& grids)
      : msubtag(subtag),
        msubsubtag(subsubtag),
        mmode(mode),
        mgrids(grids),
        mjac(target) {
    mjac.Perturbation() = perturbation;
  }
  
  const String& Subtag() const { return msubtag; }
  
  void Subtag(const String& st) { msubtag = st; }
  
  const String& SubSubtag() const { return msubsubtag; }
  
  void SubSubtag(const String& sst) { msubsubtag = sst; }
  
  const String& Mode() const { return mmode; }
  
  void Mode(const String& m) { mmode = m; }
  
  const ArrayOfVector& Grids() const { return mgrids; }
  
  void Grids(const ArrayOfVector& g) { mgrids = g; }
 
  Index nelem() const {
    Index i = 1;
    for (Index j = 0; j < mgrids.nelem(); ++j) {
      i *= mgrids[j].nelem();
    }
    return i;
  }
  
  Jacobian::Target& Target() {return mjac;}
  
  const Jacobian::Target& Target() const {return mjac;}
  
  void Target(const Jacobian::Target& jac) {mjac=jac;}
  
  const QuantumIdentifier& QuantumIdentity() const {
    return mjac.QuantumIdentity();
  }
  
  Jacobian::Line LineType() const noexcept {return mjac.LineType();}
  
  bool operator==(Jacobian::Atm other) const noexcept {return mjac==other;}
  
  bool operator==(Jacobian::Line other) const noexcept {return mjac==other;}
  
  bool operator==(Jacobian::Sensor other) const noexcept {return mjac==other;}
  
  bool operator==(Jacobian::Special other) const noexcept {return mjac==other;}
  
  bool operator==(Jacobian::Type other) const noexcept {return mjac==other;}
  
  bool operator==(const ArrayOfSpeciesTag& st) const noexcept {return mjac==Jacobian::Special::ArrayOfSpeciesTagVMR and mjac.SpeciesList() == st;}
  
  void QuantumIdentity(const QuantumIdentifier& qi) { mjac.QuantumIdentity() = qi; }
  
  bool propmattype() const noexcept {
    return mjac == Jacobian::Type::Line or 
           mjac == Jacobian::Type::Atm or 
           mjac == Jacobian::Special::ArrayOfSpeciesTagVMR;
  }
  
  bool is_wind() const noexcept {return mjac.isWind();}
  
  bool is_mag() const noexcept {return mjac.isMagnetic();}
  
  void SetTransformationFunc(const String& s) { transformation_func = s; }
  void SetTFuncParameters(const Vector& p) { tfunc_parameters = p; }
  void SetTransformationMatrix(const Matrix& A) { transformation_matrix = A; }
  void SetOffsetVector(const Vector& b) { offset_vector = b; }
  bool HasAffine() const { return !transformation_matrix.empty(); }
  const String& TransformationFunc() const { return transformation_func; }
  const Vector& TFuncParameters() const { return tfunc_parameters; }
  const Matrix& TransformationMatrix() const { return transformation_matrix; }
  const Vector& OffsetVector() const { return offset_vector; }
 
  bool HasSameInternalsAs(const RetrievalQuantity& a) const {
    return a.msubtag == msubtag and
           a.msubsubtag == msubsubtag and a.mmode == mmode and
           a.mjac.sameTargetType(mjac);
  }
  
  String& SubTag() {return msubtag;}
  String& SubSubTag() {return msubsubtag;}
  String& Mode() {return mmode;}
  ArrayOfVector& Grids() {return mgrids;}
  String& TransformationFunc() {return transformation_func;}
  Vector& TFuncParameters() {return tfunc_parameters;}
  Matrix& Transformation() {return transformation_matrix;}
  Vector& Offset() {return offset_vector;}
 
 private:
  String msubtag;
  String msubsubtag;
  String mmode;
  ArrayOfVector mgrids;
  Jacobian::Target mjac;
 
  String transformation_func;
  Vector tfunc_parameters;
 
  Matrix transformation_matrix;
  Vector offset_vector;
};
 
ostream& operator<<(ostream& os, const RetrievalQuantity& ot);
 
typedef Array<RetrievalQuantity> ArrayOfRetrievalQuantity;
 
// A macro to loop analytical jacobian quantities
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)                             \
  for (Index iq = 0; iq < jacobian_quantities.nelem(); iq++) {              \
    if (not(jacobian_quantities[iq] == Jacobian::Type::Sensor) and          \
        not(jacobian_quantities[iq] == Jacobian::Special::SurfaceString)) { \
      what_to_do                                                            \
    }                                                                       \
  }
// A macro to loop analytical jacobian quantities
#define FOR_ANALYTICAL_JACOBIANS_DO2(what_to_do)                  \
  for (Index iq = 0; iq < jacobian_quantities.nelem(); iq++) {    \
    if (not(jacobian_quantities[iq] == Jacobian::Type::Sensor)) { \
      what_to_do                                                  \
    }                                                             \
  }
 
//======================================================================
//             Index ranges and transformation functions
//======================================================================
 
void jac_ranges_indices(ArrayOfArrayOfIndex& jis,
                        bool& any_affine,
                        const ArrayOfRetrievalQuantity& jqs,
                        const bool& before_affine = false);
 
void transform_jacobian(Matrix& jacobian,
                        const Vector x,
                        const ArrayOfRetrievalQuantity& jqs);
 
void transform_x(Vector& x, const ArrayOfRetrievalQuantity& jqs);
 
void transform_x_back(Vector& x_t,
                      const ArrayOfRetrievalQuantity& jqs,
                      bool revert_functional_transforms = true);
 
//======================================================================
//             Functions related to calculation of Jacobian
//======================================================================
 
bool check_retrieval_grids(ArrayOfVector& grids,
                           ostringstream& os,
                           const Vector& p_grid,
                           const Vector& lat_grid,
                           const Vector& lon_grid,
                           const Vector& p_retr,
                           const Vector& lat_retr,
                           const Vector& lon_retr,
                           const String& p_retr_name,
                           const String& lat_retr_name,
                           const String& lon_retr_name,
                           const Index& dim);
 
 
bool check_retrieval_grids(ArrayOfVector& grids,
                           ostringstream& os,
                           const Vector& lat_grid,
                           const Vector& lon_grid,
                           const Vector& lat_retr,
                           const Vector& lon_retr,
                           const String& lat_retr_name,
                           const String& lon_retr_name,
                           const Index& dim);
 
void diy_from_path_to_rgrids(Tensor3View diy_dx,
                             const RetrievalQuantity& jacobian_quantity,
                             ConstTensor3View diy_dpath,
                             const Index& atmosphere_dim,
                             const Ppath& ppath,
                             ConstVectorView ppath_p);
 
void diy_from_pos_to_rgrids(Tensor3View diy_dx,
                            const RetrievalQuantity& jacobian_quantity,
                            ConstMatrixView diy_dpos,
                            const Index& atmosphere_dim,
                            ConstVectorView rtp_pos);
 
ArrayOfTensor3 get_standard_diy_dpath(const ArrayOfRetrievalQuantity& jacobian_quantities, Index np, Index nf, Index ns, bool active);
 
ArrayOfTensor3 get_standard_starting_diy_dx(const ArrayOfRetrievalQuantity& jacobian_quantities, Index np, Index nf, Index ns, bool active);
 
ArrayOfIndex get_pointers_for_analytical_species(const ArrayOfRetrievalQuantity& jacobian_quantities,
                                                 const ArrayOfArrayOfSpeciesTag& abs_species);
 
ArrayOfIndex get_pointers_for_scat_species(const ArrayOfRetrievalQuantity& jacobian_quantities,
                                           const ArrayOfString& scat_species,
                                           const bool cloudbox_on);
 
template <std::size_t N>
Index do_analytical_jacobian(const ArrayOfRetrievalQuantity& jacobian_quantities) {
  static_assert(N == 1 or N == 2, "FOR_ANALYTICAL_JACOBIANS_DO or FOR_ANALYTICAL_JACOBIANS_DO2");
  if constexpr (N == 1) FOR_ANALYTICAL_JACOBIANS_DO(return 1;)
  else if constexpr (N == 2) FOR_ANALYTICAL_JACOBIANS_DO2(return 1;)
  return 0;
}
 
void jacobian_type_extrapol(ArrayOfGridPos& gp);
 
void polynomial_basis_func(Vector& b, const Vector& x, const Index& poly_coeff);
 
void calcBaselineFit(Vector& y_baseline,
                     const Vector& x,
                     const Index& mblock_index,
                     const Sparse& sensor_response,
                     const ArrayOfIndex& sensor_response_pol_grid,
                     const Vector& sensor_response_f_grid,
                     const Matrix& sensor_response_dlos_grid,
                     const RetrievalQuantity& rq,
                     const Index rq_index,
                     const ArrayOfArrayOfIndex& jacobian_indices);
 
void vmrunitscf(Numeric& x,
                const String& unit,
                const Numeric& vmr,
                const Numeric& p,
                const Numeric& t);
 
void dxdvmrscf(Numeric& x,
               const String& unit,
               const Numeric& vmr,
               const Numeric& p,
               const Numeric& t);
 
//======================================================================
//             Propmat partials descriptions
//======================================================================
 
Numeric temperature_perturbation(const ArrayOfRetrievalQuantity& js) noexcept;
 
Numeric frequency_perturbation(const ArrayOfRetrievalQuantity& js) noexcept;
 
Numeric magnetic_field_perturbation(
    const ArrayOfRetrievalQuantity& js) noexcept;
 
String propmattype_string(const RetrievalQuantity& rq);
 
//======================================================================
//             Propmat partials boolean functions
//======================================================================
 
bool is_frequency_parameter(const RetrievalQuantity& t) noexcept;
 
bool is_derived_magnetic_parameter(const RetrievalQuantity& t) noexcept;
 
bool is_nlte_parameter(const RetrievalQuantity& t) noexcept;
 
bool is_pressure_broadening_G0(const RetrievalQuantity& t) noexcept;
 
bool is_pressure_broadening_D0(const RetrievalQuantity& t) noexcept;
 
bool is_pressure_broadening_G2(const RetrievalQuantity& t) noexcept;
 
bool is_pressure_broadening_D2(const RetrievalQuantity& t) noexcept;
 
bool is_pressure_broadening_FVC(const RetrievalQuantity& t) noexcept;
 
bool is_pressure_broadening_ETA(const RetrievalQuantity& t) noexcept;
 
bool is_pressure_broadening_Y(const RetrievalQuantity& t) noexcept;
 
bool is_pressure_broadening_G(const RetrievalQuantity& t) noexcept;
 
bool is_pressure_broadening_DV(const RetrievalQuantity& t) noexcept;
 
bool is_lineshape_parameter_X0(const RetrievalQuantity& t) noexcept;
 
bool is_lineshape_parameter_X1(const RetrievalQuantity& t) noexcept;
 
bool is_lineshape_parameter_X2(const RetrievalQuantity& t) noexcept;
 
bool is_lineshape_parameter_bar_linemixing(const RetrievalQuantity& t) noexcept;
 
bool is_lineshape_parameter(const RetrievalQuantity& t) noexcept;
 
bool is_line_parameter(const RetrievalQuantity& t) noexcept;
 
bool supports_CIA(const ArrayOfRetrievalQuantity& js);
 
bool supports_hitran_xsec(const ArrayOfRetrievalQuantity& js);
 
 
bool supports_continuum(const ArrayOfRetrievalQuantity& js);
 
bool supports_relaxation_matrix(const ArrayOfRetrievalQuantity& js);
 
bool supports_lookup(const ArrayOfRetrievalQuantity& js);
 
bool supports_zeeman(const ArrayOfRetrievalQuantity& js);
 
bool supports_faraday(const ArrayOfRetrievalQuantity& js);
 
bool supports_propmat_clearsky(const ArrayOfRetrievalQuantity& js);
 
bool species_match(const RetrievalQuantity& rq, const ArrayOfSpeciesTag& ast);
 
bool species_match(const RetrievalQuantity& rq, const Species::Species species);
 
bool species_iso_match(const RetrievalQuantity& rq,
                       const Species::IsotopeRecord& ir);
 
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity& js) noexcept;
 
struct jacobianVMRcheck {
  bool test;
  const QuantumIdentifier& qid;
};
 
jacobianVMRcheck do_vmr_jacobian(const ArrayOfRetrievalQuantity& js,
                                 const QuantumIdentifier& line_qid) noexcept;
 
bool do_line_center_jacobian(const ArrayOfRetrievalQuantity& js) noexcept;
 
bool do_frequency_jacobian(const ArrayOfRetrievalQuantity& js) noexcept;
 
bool do_magnetic_jacobian(const ArrayOfRetrievalQuantity& js) noexcept;
 
#endif  // jacobian_h