error_analysis.h

Go to the documentation of this file.

#ifndef ERROR_ANALYSIS_H
#define ERROR_ANALYSIS_H
#include "connor_solver.h"
#include "max_a_posteriori.h"
#include "atmosphere_oco.h"
#include "forward_model.h"
#include "observation.h"
#include "fe_disable_exception.h"

namespace FullPhysics {
/****************************************************************/
class ErrorAnalysis : public Printable<ErrorAnalysis> {
public:
  ErrorAnalysis(const boost::shared_ptr<ConnorSolver>& Solver,
                const boost::shared_ptr<AtmosphereOco>& Atm,
                const boost::shared_ptr<ForwardModel>& Fm,
        const boost::shared_ptr<Observation>& inst_meas);
  ErrorAnalysis(const boost::shared_ptr<MaxAPosteriori>& Max_a_posteriori,
                const boost::shared_ptr<AtmosphereOco>& Atm,
                const boost::shared_ptr<ForwardModel>& Fm,
        const boost::shared_ptr<Observation>& inst_meas);
  ErrorAnalysis(const boost::shared_ptr<ConnorSolver>& Solver,
                const boost::shared_ptr<RtAtmosphere>& Atm,
                const boost::shared_ptr<ForwardModel>& Fm,
        const boost::shared_ptr<Observation>& inst_meas);
  ErrorAnalysis(const boost::shared_ptr<MaxAPosteriori>& Max_a_posteriori,
                const boost::shared_ptr<RtAtmosphere>& Atm,
                const boost::shared_ptr<ForwardModel>& Fm,
        const boost::shared_ptr<Observation>& inst_meas);
  virtual ~ErrorAnalysis() {}

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  int num_channels() const { return fm->num_channels();}

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  blitz::Array<double, 1> modeled_radiance() const 
  { 
    if(residual().rows() == 0)
      return blitz::Array<double, 1>(0);
    return blitz::Array<double, 1>
      (residual() + meas->radiance_all().spectral_range().data()); 
  }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  double residual_sum_sq(int Band) const {
    FeDisableException disable_fp;
    boost::optional<blitz::Range> pr = fm->stacked_pixel_range(Band);
    if(!pr)
      return 0;
    blitz::Array<double, 1> res(residual());
    if(res.rows() == 0)
      return 0;
    return sum(res(*pr) * res(*pr));
  }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  double residual_mean_sq(int Band) const {
    FeDisableException disable_fp;
    boost::optional<blitz::Range> pr = fm->stacked_pixel_range(Band);
    if(!pr)
      return 0;
    return sqrt(residual_sum_sq(Band) / pr->length());
  }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  double reduced_chisq(int Band) const
  { 
    FeDisableException disable_fp;
    boost::optional<blitz::Range> pr = fm->stacked_pixel_range(Band);
    if(!pr)
      return 0;
    if(solver) {
      if(solver->residual().rows() == 0)
        return 0;
      return chisq_measure_norm(solver->residual()(*pr), 
                                solver->residual_covariance_diagonal()(*pr));
    } else {
      blitz::Array<double, 1> res(max_a_posteriori->uncert_weighted_model_measure_diff());
      if(!res.rows()) return 0;
      return sum(res(*pr)*res(*pr))/pr->length();
    }
  }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  double relative_residual_mean_sq(int Band) const 
  { 
    FeDisableException disable_fp;
    boost::optional<blitz::Range> pr = fm->stacked_pixel_range(Band);
    if(!pr)
      return 0;
    double result = residual_mean_sq(Band);
    return result ? (result/signal(Band)):0;
  }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  double chisq_measure_norm(const blitz::Array<double, 1>& Residual,
            const blitz::Array<double, 1>& Residual_cov_diag) const
  { 
    FeDisableException disable_fp;
    if (residual().rows() == 0)
      return 0;
    return sum(Residual * Residual / Residual_cov_diag) / Residual.rows(); 
  }

  double signal(int band) const;
  double noise(int band) const;

  double xco2_measurement_error() const; 
  double xco2_smoothing_error() const;
  double xco2_uncertainty() const;
  double xco2_interference_error() const;
  blitz::Array<double, 1> xco2_gain_vector() const;

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  double gamma_last_step() const 
  { return (solver ? solver->gamma_last_step() : -1); }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  double xco2_uncert_noise() const
  { 
      FeDisableException disable_fp;
      return sqrt(xco2_measurement_error());
  }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  double xco2_uncert_smooth() const
  { 
      FeDisableException disable_fp;
      return sqrt(xco2_smoothing_error());
  }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  double xco2_uncert_interf() const
  {
      FeDisableException disable_fp;
      return sqrt(xco2_interference_error());
  }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  double degrees_of_freedom_full_vector() const
  { 
      FeDisableException disable_fp;
      return sum(averaging_kernel()(i1, i1));
  }

//-----------------------------------------------------------------------
//-----------------------------------------------------------------------

  double degrees_of_freedom_xco2() const 
  { 
      FeDisableException disable_fp;
      return sum(where(xco2_state_used(),
                     averaging_kernel()(i1, i1), 0));
  }

  blitz::Array<double, 1> xco2_avg_kernel() const;
  blitz::Array<double, 2> co2_averaging_kernel() const;
  blitz::Array<double, 1> xco2_avg_kernel_full() const;
  blitz::Array<double, 1> xco2_avg_kernel_norm() const;
  blitz::Array<double, 1> xco2_correlation_interf() const;
  blitz::Array<double, 1> interference_smoothing_uncertainty() const;

  void print(std::ostream& Os) const { Os << "ErrorAnalysis";}
private:
  blitz::Array<double, 1> residual() const
  {
    if(solver)
      return solver->residual();
    else
      return max_a_posteriori->model_measure_diff();
  }
  blitz::Array<double, 2> averaging_kernel() const
  {
    if(solver)
      return solver->averaging_kernel();
    else
      return max_a_posteriori->averaging_kernel();
  }
  blitz::Array<double, 2> aposteriori_covariance() const
  {
    if(solver)
      return solver->aposteriori_covariance();
    else
      return max_a_posteriori->a_posteriori_covariance();
  }
  blitz::Array<double, 2> apriori_covariance() const
  {
    if(solver)
      return solver->apriori_covariance();
    else
      return max_a_posteriori->a_priori_cov();
  }
  blitz::Array<bool, 1> xco2_state_used() const 
  { return atm->absorber_ptr()->absorber_vmr("CO2")->state_used(); }
  AutoDerivative<double> xco2() const 
  { return atm->absorber_ptr()->xgas("CO2"); }
  blitz::Array<double, 1> dxco2_dstate() const;
  // Only one of solver or max_a_posteriori will be nonnull.
  boost::shared_ptr<ConnorSolver> solver;
  boost::shared_ptr<MaxAPosteriori> max_a_posteriori;
  boost::shared_ptr<AtmosphereOco> atm;
  boost::shared_ptr<ForwardModel> fm;
  boost::shared_ptr<Observation> meas;
  blitz::Array<double, 2> hmat() const;
  blitz::Array<double, 2> ht_c_h() const;
  // Used in a lot of places, so define once here.
  blitz::firstIndex i1; blitz::secondIndex i2; blitz::thirdIndex i3; 
  blitz::fourthIndex i4;
};
}
#endif