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 <iostream>
#include <map>
#include <stdexcept>
#include "abs_species_tags.h"
#include "agenda_class.h"
#include "array.h"
#include "bifstream.h"
#include "interpolation.h"
#include "logic.h"
#include "matpackI.h"
#include "methods.h"
#include "mystring.h"
#include "ppath.h"

#include "quantum.h"

enum class JacPropMatType : Index {
  VMR,
  Electrons,
  Particulates,
  Temperature,

  MagneticMagnitude,
  MagneticU,
  MagneticV,
  MagneticW,

  WindMagnitude,
  WindU,
  WindV,
  WindW,
  Frequency,  // Note:  This is how wind is done internal to propmat agenda

  LineStrength,
  LineCenter,
  
  LineSpecialParameter1,

  LineShapeG0X0,
  LineShapeG0X1,
  LineShapeG0X2,
  LineShapeG0X3,

  LineShapeD0X0,
  LineShapeD0X1,
  LineShapeD0X2,
  LineShapeD0X3,

  LineShapeG2X0,
  LineShapeG2X1,
  LineShapeG2X2,
  LineShapeG2X3,

  LineShapeD2X0,
  LineShapeD2X1,
  LineShapeD2X2,
  LineShapeD2X3,

  LineShapeFVCX0,
  LineShapeFVCX1,
  LineShapeFVCX2,
  LineShapeFVCX3,

  LineShapeETAX0,
  LineShapeETAX1,
  LineShapeETAX2,
  LineShapeETAX3,

  LineShapeYX0,
  LineShapeYX1,
  LineShapeYX2,
  LineShapeYX3,

  LineShapeGX0,
  LineShapeGX1,
  LineShapeGX2,
  LineShapeGX3,

  LineShapeDVX0,
  LineShapeDVX1,
  LineShapeDVX2,
  LineShapeDVX3,

  NLTE,

  NotPropagationMatrixType
};


class RetrievalQuantity {
 public:
  RetrievalQuantity()
      : mmaintag(),
        msubtag(),
        msubsubtag(),
        mmode(),
        manalytical(-1),
        mperturbation(0.),
        mgrids(),
        mquantumidentifier(),
        mproptype(JacPropMatType::NotPropagationMatrixType),
        mintegration_flag(false) { /* Nothing to do here. */
  }

  RetrievalQuantity(const String& maintag,
                    const String& subtag,
                    const String& subsubtag,
                    const String& mode,
                    const Index& analytical,
                    const Numeric& perturbation,
                    const ArrayOfVector& grids)
      : mmaintag(maintag),
        msubtag(subtag),
        msubsubtag(subsubtag),
        mmode(mode),
        manalytical(analytical),
        mperturbation(perturbation),
        mgrids(grids),
        mquantumidentifier(),
        mproptype(JacPropMatType::NotPropagationMatrixType),
        mintegration_flag(false) {
    // With Matpack, initialization of mgrids from grids should work correctly.
  }

  const String& MainTag() const { return mmaintag; }
  
  void MainTag(const String& mt) { mmaintag = mt; }
  
  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 Index& Analytical() const { return manalytical; }
  
  void Analytical(const Index& m) { manalytical = m; }
  
  const Numeric& Perturbation() const { return mperturbation; }
  
  void Perturbation(const Numeric& p) { mperturbation = p; }
  
  const ArrayOfVector& Grids() const { return mgrids; }
  
  void Grids(const ArrayOfVector& g) { mgrids = g; }
  
  JacPropMatType PropMatType() const { return mproptype; }
  
  void PropType(const JacPropMatType t) { mproptype = t; }
  
  bool operator==(const JacPropMatType t) const { return t == mproptype; }
  
  bool operator!=(const JacPropMatType t) const { return not operator==(t); }

  Index nelem() const {
    Index i = 1;
    for (Index j = 0; j < mgrids.nelem(); ++j) {
      i *= mgrids[j].nelem();
    }
    return i;
  }
  
  const QuantumIdentifier& QuantumIdentity() const {
    return mquantumidentifier;
  }
  
  void QuantumIdentity(const QuantumIdentifier& qi) { mquantumidentifier = qi; }

  bool Integration() const { return mintegration_flag; }
  
  void IntegrationOn() { mintegration_flag = true; }
  
  void IntegrationOff() { mintegration_flag = false; }

  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.mmaintag == mmaintag and a.msubtag == msubtag and
           a.msubsubtag == msubsubtag and a.mmode == mmode and
           a.manalytical == manalytical and
           a.mquantumidentifier == mquantumidentifier and
           a.mproptype == mproptype;
  }
  
  String& MainTag() {return mmaintag;}
  String& SubTag() {return msubtag;}
  String& SubSubTag() {return msubsubtag;}
  String& Mode() {return mmode;}
  Index& Analytical() {return manalytical;}
  Numeric& Perturbation() {return mperturbation;}
  ArrayOfVector& Grids() {return mgrids;}
  QuantumIdentifier& QuantumIdentity() {return mquantumidentifier;}
  JacPropMatType Proptype() const {return mproptype;}
  Index Proptype(JacPropMatType x) {if (validProptype(x)) {mproptype = x; return EXIT_SUCCESS;} else return EXIT_FAILURE;}
  static bool validProptype(JacPropMatType x) noexcept {return Index(x) <= Index(JacPropMatType::NotPropagationMatrixType) and Index(x) >= 0;}
  void Integration(bool x) {mintegration_flag = x;}
  String& TransformationFunc() {return transformation_func;}
  Vector& TFuncParameters() {return tfunc_parameters;}
  Matrix& Transformation() {return transformation_matrix;}
  Vector& Offset() {return offset_vector;}

 private:
  String mmaintag;
  String msubtag;
  String msubsubtag;
  String mmode;
  Index manalytical;
  Numeric mperturbation;
  ArrayOfVector mgrids;
  QuantumIdentifier mquantumidentifier;
  JacPropMatType mproptype;
  bool mintegration_flag;

  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 (jacobian_quantities[iq].Analytical()) {                \
      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 (jacobian_quantities[iq].Analytical() ||                 \
        jacobian_quantities[iq].MainTag() == SURFACE_MAINTAG) { \
      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);
 
void get_pointers_for_analytical_jacobians(
    ArrayOfIndex& abs_species_i,
    ArrayOfIndex& scat_species_i,
    ArrayOfIndex& is_t,
    ArrayOfIndex& wind_i,
    ArrayOfIndex& magfield_i,
    const ArrayOfRetrievalQuantity& jacobian_quantities,
    const ArrayOfArrayOfSpeciesTag& abs_species,
    const Index& cloudbox_on,
    const ArrayOfString& scat_species);

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);

enum class JacobianType : Index {
  None = 0,  // Setting to nil means that (bool)0 and (bool)N still works.
  Temperature,
  WindFieldU,
  WindFieldV,
  WindFieldW,
  AbsWind,
  MagFieldU,
  MagFieldV,
  MagFieldW,
  AbsMag,
  Other
};

void get_diydx(VectorView diy1,
               VectorView diy2,
               ConstMatrixView ImT,
               ConstMatrixView cumulative_transmission,
               ConstMatrixView dT1,
               ConstMatrixView dT2,
               ConstVectorView iYmJ,
               ConstVectorView dJ1,
               ConstVectorView dJ2,
               const Index stokes_dim,
               const bool transmission_only = false);

//======================================================================
//             Propmat partials descriptions
//======================================================================

ArrayOfIndex equivalent_propmattype_indexes(const ArrayOfRetrievalQuantity& js);

Index equivalent_propmattype_index(const ArrayOfRetrievalQuantity& js,
                                  const Index i) noexcept;

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_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_LBL_without_phase(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_particles(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 Index species);

bool species_iso_match(const RetrievalQuantity& rq,
                       const Index species,
                       const Index iso);

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