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 "arts_conversions.h"
#include "arts_options.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_struct.h"
#include "quantum_numbers.h"  
#include "species_tags.h"
#include <iostream>
#include <map>
#include <stdexcept>
#include <utility>
 
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,
          ECS_SCALINGX0, ECS_SCALINGX1, ECS_SCALINGX2, ECS_SCALINGX3,
          ECS_BETAX0, ECS_BETAX1, ECS_BETAX2, ECS_BETAX3,
          ECS_LAMBDAX0, ECS_LAMBDAX1, ECS_LAMBDAX2, ECS_LAMBDAX3,
          ECS_DCX0, ECS_DCX1, ECS_DCX2, ECS_DCX3,
          NLTE
          )
static_assert(
    Index(Line::FINAL) == 4 + 13 * Index(Options::LineShapeCoeff::FINAL),
    "Either you have added some \"Line\" parameter(s) or changed the temperature model");
 
ENUMCLASS(Sensor, char,
          FrequencyShift,
          FrequencyStretch,
          Polyfit,
          Sinefit,
          PointingZenithInterp,
          PointingZenithRecalc
          )
 
ENUMCLASS(Special, char,
          ArrayOfSpeciesTagVMR,
          ScatteringString,
          SurfaceString
          )
 
struct Target {
  Type type{Type::FINAL};
 
  Atm atm{Atm::FINAL};
 
  Line line{Line::FINAL};
 
  Sensor sensor{Sensor::FINAL};
 
  Special special{Special::FINAL};
 
  Numeric perturbation{std::numeric_limits<Numeric>::quiet_NaN()};
 
  QuantumIdentifier qid{};
 
  ArrayOfSpeciesTag species_array_id{0};
 
  String string_id{};
 
  Species::Species species_id{Species::Species::FINAL};
 
  explicit Target(Atm atype) : type(Type::Atm), atm(atype) {
    ARTS_ASSERT(good_enum(atype))
  }
 
  explicit Target(Line ltype,
                  const QuantumIdentifier& iqid,
                  Species::Species specid)
      : type(Type::Line), line(ltype), qid(iqid), species_id(specid) {
    ARTS_ASSERT(good_enum(ltype))
  }
 
  explicit Target(Sensor stype) : type(Type::Sensor), sensor(stype) {
    ARTS_ASSERT(good_enum(stype))
  }
 
  explicit Target(Special stype, ArrayOfSpeciesTag aostid)
      : type(Type::Special), special(stype), species_array_id(std::move(aostid)) {
    ARTS_ASSERT(stype == Special::ArrayOfSpeciesTagVMR,
                "Only for Special::ArrayOfSpeciesTagVMR, but you fed: ",
                special)
  }
 
  explicit Target(Special stype, String sid)
      : type(Type::Special), special(stype), string_id(std::move(sid)) {
    ARTS_ASSERT(
        stype == Special::SurfaceString or stype == Special::ScatteringString,
        "Only for Special::ScatteringString or Special::SurfaceString, but you fed: ",
        special)
  }
 
  explicit Target() = default;
 
  bool operator==(Atm other) const noexcept { return other == atm; }
 
  bool operator==(Line other) const noexcept { return other == line; }
 
  bool operator==(Sensor other) const noexcept { return other == sensor; }
 
  bool operator==(Special other) const noexcept { return other == special; }
 
  bool operator==(Type other) const noexcept { return other == type; }
 
  [[nodiscard]] bool sameTargetType(const Target& other) const noexcept {
    return type == other.type and atm == other.atm and
           line == other.line and sensor == other.sensor and
           special == other.special;
  }
 
  [[nodiscard]] std::string_view TargetType() const noexcept {
    return toString(type);
  }
 
  void TargetType(const std::string_view& s) noexcept { type = toType(s); }
 
  void TargetSubType(const std::string_view& s) noexcept {
    atm = Atm::FINAL;
    line = Line::FINAL;
    sensor = Sensor::FINAL;
    special = Special::FINAL;
 
    switch (type) {
      case Type::Atm:
        atm = toAtm(s);
        break;
      case Type::Line:
        line = toLine(s);
        break;
      case Type::Sensor:
        sensor = toSensor(s);
        break;
      case Type::Special:
        special = toSpecial(s);
        break;
      case Type::FINAL: { /* leave last, don't use default */
      }
    }
  }
 
  [[nodiscard]] std::string_view TargetSubType() const noexcept {
    switch (type) {
      case Type::Atm:
        return toString(atm);
      case Type::Line:
        return toString(line);
      case Type::Sensor:
        return toString(sensor);
      case Type::Special:
        return toString(special);
      case Type::FINAL: { /* leave last, don't use default */
      }
    }
    return "BAD SUBTYPE";
  }
 
  [[nodiscard]] bool TargetTypeOK() const noexcept { return good_enum(type); }
 
  [[nodiscard]] 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(special) + good_enum(sensor) + good_enum(line) +
              good_enum(atm))) {
      switch (type) {
        case Type::Special:
          return good_enum(special);
        case Type::Sensor:
          return good_enum(sensor);
        case Type::Line:
          return good_enum(line);
        case Type::Atm:
          return good_enum(atm);
        case Type::FINAL: { /* leave last, don't use default */
        }
      }
    }
 
    return false;
  }
 
  [[nodiscard]] bool isSpeciesVMR() const noexcept {
    return line == Line::VMR or special == Special::ArrayOfSpeciesTagVMR;
  }
 
  [[nodiscard]] bool isWind() const noexcept {
    return atm == Atm::WindMagnitude or atm == Atm::WindU or
           atm == Atm::WindV or atm == Atm::WindW;
  }
 
  [[nodiscard]] bool isMagnetic() const noexcept {
    return atm == Atm::MagneticMagnitude or atm == Atm::MagneticU or
           atm == Atm::MagneticV or atm == Atm::MagneticW;
  }
 
  [[nodiscard]] bool isFrequency() const noexcept {
    return sensor == Sensor::FrequencyStretch or
           sensor == Sensor::FrequencyShift;
  }
 
  [[nodiscard]] bool isPointing() const noexcept {
    return sensor == Sensor::PointingZenithInterp or
           sensor == Sensor::PointingZenithRecalc;
  }
 
  [[nodiscard]] bool needQuantumIdentity() const noexcept {
    return type == Type::Line;
  }
 
  [[nodiscard]] bool needArrayOfSpeciesTag() const noexcept {
    return special == Special::ArrayOfSpeciesTagVMR;
  }
 
  [[nodiscard]] bool needString() const noexcept {
    return special == Special::ScatteringString or
           special == Special::SurfaceString;
  }
 
  friend std::ostream& operator<<(std::ostream& os, const Target& x);
};  // Target
} // namespace Jacobian
using JacobianTarget = Jacobian::Target;
using ArrayOfJacobianTarget = Array<Jacobian::Target>;
 
 
class RetrievalQuantity {
 public:
  RetrievalQuantity()
      : msubtag(),
        msubsubtag(),
        mmode(),
        mgrids(),
        mjac() { /* Nothing to do here. */
  }
 
  RetrievalQuantity(Jacobian::Target  target,
                    String  subtag,
                    String  subsubtag,
                    String  mode,
                    const Numeric& perturbation,
                    ArrayOfVector  grids)
      : msubtag(std::move(subtag)),
        msubsubtag(std::move(subsubtag)),
        mmode(std::move(mode)),
        mgrids(std::move(grids)),
        mjac(std::move(target)) {
    mjac.perturbation = perturbation;
  }
  
  [[nodiscard]] const String& Subtag() const { return msubtag; }
  
  void Subtag(const String& st) { msubtag = st; }
  
  [[nodiscard]] const String& SubSubtag() const { return msubsubtag; }
  
  void SubSubtag(const String& sst) { msubsubtag = sst; }
  
  [[nodiscard]] const String& Mode() const { return mmode; }
  
  void Mode(const String& m) { mmode = m; }
  
  [[nodiscard]] const ArrayOfVector& Grids() const { return mgrids; }
  
  void Grids(const ArrayOfVector& g) { mgrids = g; }
 
  [[nodiscard]] 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;}
  
  [[nodiscard]] const Jacobian::Target& Target() const {return mjac;}
  
  void Target(const Jacobian::Target& jac) {mjac=jac;}
  
  [[nodiscard]] const QuantumIdentifier& QuantumIdentity() const {
    return mjac.qid;
  }
  
  [[nodiscard]] Jacobian::Line LineType() const noexcept {return mjac.line;}
  
  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.species_array_id == st;}
  
  void QuantumIdentity(const QuantumIdentifier& qi) { mjac.qid = qi; }
  
  [[nodiscard]] bool propmattype() const noexcept {
    return mjac == Jacobian::Type::Line or 
           mjac == Jacobian::Type::Atm or 
           mjac == Jacobian::Special::ArrayOfSpeciesTagVMR;
  }
  
  [[nodiscard]] bool is_wind() const noexcept {return mjac.isWind();}
  
  [[nodiscard]] 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; }
  [[nodiscard]] bool HasAffine() const { return !transformation_matrix.empty(); }
  [[nodiscard]] const String& TransformationFunc() const { return transformation_func; }
  [[nodiscard]] const Vector& TFuncParameters() const { return tfunc_parameters; }
  [[nodiscard]] const Matrix& TransformationMatrix() const { return transformation_matrix; }
  [[nodiscard]] const Vector& OffsetVector() const { return offset_vector; }
 
  [[nodiscard]] 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;}
 
  [[nodiscard]] const String& SubTag() const { return msubtag; }
  [[nodiscard]] const String& SubSubTag() const { return msubsubtag; }
  [[nodiscard]] const Matrix& Transformation() const { return transformation_matrix; }
  [[nodiscard]] const Vector& Offset() const { return offset_vector; }
 
  friend ostream& operator<<(ostream& os, const RetrievalQuantity& ot);
 
 private:
  String msubtag;
  String msubsubtag;
  String mmode;
  ArrayOfVector mgrids;
  Jacobian::Target mjac;
 
  String transformation_func;
  Vector tfunc_parameters;
 
  Matrix transformation_matrix;
  Vector offset_vector;
};
 
using ArrayOfRetrievalQuantity = Array<RetrievalQuantity>;
 
// 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_wind_parameter(const RetrievalQuantity& t) noexcept;
 
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_wind_jacobian(const ArrayOfRetrievalQuantity& js) noexcept;
 
bool do_frequency_jacobian(const ArrayOfRetrievalQuantity& js) noexcept;
 
bool do_magnetic_jacobian(const ArrayOfRetrievalQuantity& js) noexcept;
 
#endif  // jacobian_h