TimeDomainGIDDecon.h

#ifndef __TIME_DOMAIN_GID_DECON__
#define __TIME_DOMAIN_GID_DECON__
#include "mspass/algorithms/TimeWindow.h"
#include "mspass/algorithms/deconvolution/ComplexArray.h"
#include "mspass/algorithms/deconvolution/CNRDeconEngine.h"
#include "mspass/algorithms/deconvolution/FFTDeconOperator.h"
#include "mspass/algorithms/deconvolution/ScalarDecon.h"
#include "mspass/algorithms/deconvolution/ThreeCSpike.h"
#include "mspass/seismic/CoreTimeSeries.h"
#include "mspass/seismic/PowerSpectrum.h"
#include "mspass/seismic/Seismogram.h"
#include "mspass/seismic/TimeSeries.h"
#include "mspass/utility/AntelopePf.h"
#include "mspass/utility/Metadata.h"
#include "mspass/utility/dmatrix.h"
#include <list>
#include <memory>
#include <ostream>
#include <string>
#include <vector>
namespace mspass::algorithms::deconvolution {
class TimeDomainGIDDecon : public ScalarDecon {
public:
  TimeDomainGIDDecon(const mspass::utility::AntelopePf &md);
  TimeDomainGIDDecon(const TimeDomainGIDDecon &parent) = delete;
  TimeDomainGIDDecon &operator=(const TimeDomainGIDDecon &parent) = delete;
  void changeparameter(const mspass::utility::Metadata &md);
  int load(const mspass::seismic::CoreSeismogram &d,
           mspass::algorithms::TimeWindow dwin);
  int loadnoise(const mspass::seismic::CoreSeismogram &d,
                mspass::algorithms::TimeWindow nwin);
  int loadwavelet(const mspass::seismic::TimeSeries &wavelet);
  int loadwavelet(const mspass::seismic::CoreTimeSeries &wavelet);
  int loadnoise(const mspass::seismic::TimeSeries &noise);
  int loadnoise(const mspass::seismic::CoreTimeSeries &noise);
  int loadnoise(const mspass::seismic::PowerSpectrum &noise_spectrum);
  void clear_external_wavelet();
  void clear_external_noise();
  double deconvolution_window_start() const { return this->fftwin.start; };
  double deconvolution_window_end() const { return this->fftwin.end; };
  double wavelet_window_start() const { return this->waveletwin.start; };
  double wavelet_window_end() const { return this->waveletwin.end; };
  double output_window_start() const { return this->outputwin.start; };
  double output_window_end() const { return this->outputwin.end; };
  double noise_window_start() const { return this->nwin.start; };
  double noise_window_end() const { return this->nwin.end; };
  std::string configuration_pf_text() const { return this->config_pf_text; };
  bool leaf_parameters_have_changed() const {
    return this->leaf_parameters_changed;
  };
  mspass::utility::Metadata changed_leaf_parameters() const {
    return this->changed_leaf_metadata;
  };
  bool external_wavelet_is_loaded() const {
    return this->external_wavelet_loaded;
  };
  bool external_noise_is_loaded() const { return this->external_noise_loaded; };
  bool external_noise_spectrum_is_loaded() const {
    return this->external_noise_spectrum_loaded;
  };
  mspass::seismic::TimeSeries loaded_external_wavelet() const {
    if (!this->external_wavelet_loaded)
      return mspass::seismic::TimeSeries();
    return this->external_wavelet;
  };
  mspass::seismic::TimeSeries loaded_external_noise() const {
    if (!this->external_noise_loaded)
      return mspass::seismic::TimeSeries();
    return this->external_noise;
  };
  mspass::seismic::PowerSpectrum loaded_external_noise_spectrum() const {
    if (!this->external_noise_spectrum_loaded)
      return mspass::seismic::PowerSpectrum();
    return this->external_noise_spectrum;
  };
  int load(const mspass::seismic::CoreSeismogram &d,
           mspass::algorithms::TimeWindow dwin,
           mspass::algorithms::TimeWindow nwin);
  void process();
  ~TimeDomainGIDDecon();
  mspass::seismic::CoreSeismogram getresult();
  mspass::seismic::CoreSeismogram sparse_output();
  std::vector<double> lag_weight_vector() const;
  mspass::seismic::TimeSeries ideal_output();
  mspass::seismic::TimeSeries actual_output();
  mspass::seismic::CoreTimeSeries inverse_wavelet();
  mspass::seismic::CoreTimeSeries inverse_wavelet(double t0parent);
  mspass::utility::Metadata QCMetrics();
 
private:
  /* These are data at different stages of process.  d_all is the
  largest signal window that is assumed to have been initialized by the
  load method for this object.  d_decon is the
  result of applying the preprocessor (signal processing) deconvolution to
  all three components. d_decon is computed from d, which is a windowed
  version of the input data received by the load method.  It has a
  time duration less than or at least equal to that of d_all.
  r is the residual, which is accumulated during the iterative method.
  The time duration of r is the same as d.  It is initialized by
  convolving the inverse filter with d_all.
  n is the noise data.  It should normally be at least as long a d_all*/
  mspass::seismic::CoreSeismogram d_all, d_decon, r, n;
  /* We save the set of data lengths for clarity and a minor bit efficiency.
  ndwin is the number of samples in d_all and r.
  nnwin is the number of samples in n. */
  int ndwin, nnwin;
  /* Save the TimeWindow objects that define the extent of d_all, d_decon,
  and n.   Some things need at least some of these downstream */
  mspass::algorithms::TimeWindow dwin, outputwin, nwin, fftwin, waveletwin;
  int inverse_operator_nfft;
  int gid_noise_samples_loaded, gid_noise_samples_used;
  bool gid_noise_truncated;
  std::string config_pf_text;
  double target_dt;
  int noise_component;
  std::list<ThreeCSpike> spikes;
  std::vector<double> lag_weights;
  /* This vector contains the function time shifted and added to lag_weights
  vector after each iteration.   */
  std::vector<double> wtf;
  int nwtf; // size of wtf - cached because wtf is inside the deepest loop
 
  /* This is a pointer to the BasicDeconOperator class used for preprocessing
  Classic use of inheritance to simplify the api. */
  std::unique_ptr<ScalarDecon> preprocessor;
  std::unique_ptr<CNRDeconEngine> cnrprocessor;
  mspass::seismic::TimeSeries current_wavelet;
  mspass::seismic::TimeSeries external_wavelet;
  mspass::seismic::TimeSeries external_noise;
  mspass::seismic::PowerSpectrum external_noise_spectrum;
  std::vector<std::vector<double>> ns_noise_components;
  bool external_wavelet_loaded, external_noise_loaded,
      external_noise_spectrum_loaded, external_wavelet_allowed;
  bool processed;
  bool residual_noise_from_external;
  bool leaf_parameters_changed;
  int gid_analysis_samples, gid_wavelet_samples, gid_alignment_offset_samples;
  double gid_analysis_t0, gid_wavelet_t0;
  mspass::utility::Metadata changed_leaf_metadata;
  mspass::utility::Metadata leaf_operator_metadata;
 
  /* Scalar leaves expose a normalized actual_output, while their process
   * result retains the physical inverse gain.  These values record the raw
   * source-through-inverse zero-lag gain and the reciprocal applied to every
   * inverse-domain quantity used by GID. */
  double gid_leaf_raw_zero_lag_gain;
  double gid_inverse_domain_amplitude_scale;
 
  /* This parameter is set in the constructor.  It would normally be half the
  length of the fir representation of the inverse wavelet.*/
  int wavelet_pad;
  std::vector<double> actual_o_fir;
  int actual_o_0; // offset from sample zero for zero lag position
  IterDeconType decon_type;
  /* This is called by the constructor to create the wtf penalty function */
  void construct_weight_penalty_function(const mspass::utility::Metadata &md);
  void invalidate_processing_state();
  void ensure_inverse_operator_size(const int data_npts,
                                    const int wavelet_npts,
                                    const int noise_npts);
  int actual_inverse_operator_size() const;
  void update_residual_matrix(ThreeCSpike spk);
  void update_lag_weights(int col, const double candidate_amplitude);
  double compute_resid_linf_floor(
      const mspass::seismic::CoreSeismogram &noise);
  bool has_not_converged();
  /* These are convergence attributes.   lw_inf indicates Linf norm of
  lag_weight array, lw_l2 is L2 metric of lag_weight, resid_inf is Linf
  norm of residual vector, and resid_l2 is L2 of resid matrix.   prev
  modifier means the size of that quantity in the previous iteration.  initial
  means initial value at the top of the loop.*/
  double lw_linf_initial, lw_linf_prev;
  double lw_l2_initial, lw_l2_prev;
  double resid_linf_initial, resid_linf_prev;
  double resid_l2_initial, resid_l2_prev;
  /* These are convergence parameters for the different tests */
  int iter_count, iter_max; // actual iteration count and ceiling to break loop
  /*lw metrics are scaled with range of 0 to 1.  l2 gets scaled by number of
  points and so can use a similar absolute scale. */
  double lw_linf_floor, lw_l2_floor;
  /* We use a probability level to define the floor Linf of the residual
  matrix.   For L2 we use the conventional fractional improvement metric. */
  double resid_linf_prob, resid_linf_floor;
  double resid_l2_tol;
  double ns_peak_sigma_threshold, ns_peak_probability_threshold;
  double ns_residual_noise_ratio_floor, ns_peak_threshold;
  double ns_last_peak_significance, ns_noise_l2, ns_noise_amplitude_rms;
  double ns_noise_component_sigma_rms, ns_noise_component_sigma_rms_robust;
  double ns_noise_component_rms_aggregate;
  bool ns_noise_component_sigma_rms_fallback_used;
  double ns_residual_rms_initial, ns_residual_rms_final;
  double ns_peak_threshold_empirical, ns_peak_threshold_sigma;
  double ns_noise_amplitude_robust, ns_last_candidate_amplitude;
  int ns_noise_samples_at_or_above_peak_threshold;
  int ns_noise_amplitude_sample_count;
  int ns_initial_stationary_null_search_lag_count;
  double ns_initial_stationary_null_expected_noise_exceedances;
  int ns_last_selected_candidate_lag;
  double ns_last_selected_candidate_lag_weight;
  double ns_last_selected_candidate_weighted_amplitude;
  double ns_max_raw_candidate_amplitude, ns_max_raw_candidate_significance;
  int ns_max_raw_candidate_lag;
  bool ns_last_scan_raw_significant_candidate_remaining;
  double ns_final_scan_max_raw_candidate_amplitude,
      ns_final_scan_max_raw_candidate_significance;
  int ns_final_scan_max_raw_candidate_lag;
  bool ns_final_scan_raw_significant_candidate_remaining;
  double ns_final_scan_existing_support_max_raw_amplitude,
      ns_final_scan_existing_support_max_raw_significance;
  int ns_final_scan_existing_support_max_raw_lag;
  int ns_final_scan_significant_candidate_count;
  int ns_final_scan_best_trial_lag;
  double ns_final_scan_best_trial_residual_l2,
      ns_final_scan_best_trial_fractional_improvement;
  int ns_final_scan_decision_candidate_lag,
      ns_final_scan_global_acceptable_candidate_count;
  double ns_final_scan_decision_trial_residual_l2,
      ns_final_scan_decision_trial_fractional_improvement;
  std::string ns_final_scan_decision;
  bool ns_final_scan_acceptable_candidate_remaining;
  std::vector<double> ns_noise_component_rms;
  std::vector<int> ns_candidate_lag_history, ns_candidate_accepted_history;
  std::vector<double> ns_candidate_lag_time_history,
      ns_candidate_amplitude_history,
      ns_candidate_threshold_history, ns_candidate_significance_history;
  std::vector<double> ns_candidate_post_residual_rms_ratio_history;
  std::vector<double> ns_candidate_residual_l2_before_history,
      ns_candidate_trial_residual_l2_history,
      ns_candidate_post_refit_residual_l2_history,
      ns_candidate_fractional_improvement_history,
      ns_candidate_state_fractional_improvement_history;
  std::vector<int> ns_candidate_periodic_refit_applied_history,
      ns_candidate_final_refit_applied_history,
      ns_candidate_trial_evaluated_history,
      ns_candidate_metric_available_history;
  std::vector<std::string> ns_candidate_stop_history;
  /* Compatibility counters for legacy pre-trial candidate scans. */
  int legacy_eq15_candidates_tested, legacy_eq15_candidates_rejected;
  /* Counters for the actual post-acceptance Eq.(15) state comparison. */
  int legacy_eq15_post_acceptance_state_tests,
      legacy_eq15_post_acceptance_floor_stops;
  int legacy_eq15_candidates_below_floor, legacy_eq15_candidates_non_decreasing,
      legacy_eq15_candidates_nonfinite,
      legacy_eq15_rejected_lag_samples_truncated,
      legacy_eq15_rejected_iteration_samples_truncated;
  double legacy_eq15_last_trial_fractional_improvement;
  std::string legacy_eq15_stop_detail;
  std::vector<double> legacy_eq15_rejected_lag_times;
  std::vector<int> legacy_eq15_rejected_candidates_per_iteration;
  int ns_max_spikes, ns_refit_interval;
  /* Audit post-refit terminal scans and any continuation they unlock. */
  int ns_refit_epochs, ns_refit_resume_count;
  double ns_ridge_beta, ns_fractional_improvement_final,
      ns_fractional_improvement_state_final;
  bool ns_use_empirical_noise_threshold, ns_converged;
  bool ns_final_refit_applied;
  std::string ns_stop_reason, ns_provisional_stop_reason_before_final_refit;
  bool gid_converged;
  std::string gid_stop_reason;
  std::string lag_weight_penalty_function;
  double lag_weight_penalty_scale_factor;
  int lag_weight_function_width;
  std::vector<double> adaptive_penalty_memory;
  std::vector<double> adaptive_penalty_retention;
  double adaptive_penalty_last_confidence;
  double adaptive_penalty_last_immediate_strength;
  double adaptive_penalty_last_specificity;
  double adaptive_penalty_last_decay_factor;
  double adaptive_penalty_noise_amplitude;
  double adaptive_penalty_memory_linf;
  double adaptive_penalty_memory_l2;
  double group_sparse_lambda, group_sparse_lambda_scale;
  double group_sparse_lambda_used, group_sparse_tolerance;
  double group_sparse_active_threshold, group_sparse_active_threshold_scale;
  double group_sparse_active_threshold_quantile;
  double group_sparse_active_threshold_quantile_value;
  double group_sparse_active_threshold_used;
  double group_sparse_objective_initial, group_sparse_objective_final;
  double group_sparse_fractional_improvement_final;
  double group_sparse_debiased_objective_final;
  double group_sparse_debiased_fractional_improvement_final;
  double group_sparse_refit_gram_condition_number;
  double group_sparse_refit_relative_ridge_beta;
  double group_sparse_refit_residual_l2_pre, group_sparse_refit_residual_l2_post;
  double group_sparse_refit_maximum_amplitude_pre;
  double group_sparse_refit_maximum_amplitude_post;
  int group_sparse_max_iterations, group_sparse_iterations;
  int group_sparse_active_groups;
  bool group_sparse_converged, group_sparse_refit_condition_guard_applied;
  bool group_sparse_refit_fallback_to_pre_debias;
  std::string group_sparse_refit_fallback_reason;
 
};
} // namespace mspass::algorithms::deconvolution
#endif