FrequencyDomainGIDDecon.h

#ifndef __FREQUENCY_DOMAIN_GID_DECON__
#define __FREQUENCY_DOMAIN_GID_DECON__
#include "mspass/algorithms/TimeWindow.h"
#include "mspass/algorithms/deconvolution/CNRDeconEngine.h"
#include "mspass/algorithms/deconvolution/FFTDeconOperator.h"
#include "mspass/algorithms/deconvolution/ScalarDecon.h"
#include "mspass/algorithms/deconvolution/ShapingWavelet.h"
#include "mspass/algorithms/deconvolution/ThreeCSpike.h"
#include "mspass/seismic/CoreSeismogram.h"
#include "mspass/seismic/CoreTimeSeries.h"
#include "mspass/seismic/PowerSpectrum.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 <string>
#include <vector>
 
namespace mspass::algorithms::deconvolution {
class FrequencyDomainGIDDecon : public ScalarDecon {
public:
  FrequencyDomainGIDDecon(const mspass::utility::AntelopePf &md);
  FrequencyDomainGIDDecon(const FrequencyDomainGIDDecon &parent) = delete;
  FrequencyDomainGIDDecon &operator=(const FrequencyDomainGIDDecon &parent) =
      delete;
  ~FrequencyDomainGIDDecon();
  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();
  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:
  mspass::seismic::CoreSeismogram d_all, d_decon, r, n;
  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 ndwin, nnwin, noise_component;
  int actual_o_0, iter_count, iter_max;
  double residual_ratio_floor, residual_improvement_floor;
  double resid_l2_initial, resid_l2_prev, resid_l2_final;
  double resid_linf_initial, resid_linf_final;
  double lag_weight_linf_final, lag_weight_l2_final;
  IterDeconType decon_type;
  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;
  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;
  std::vector<double> actual_o_fir;
  /* 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;
  std::vector<double> lag_weights, lag_weight_penalty;
  std::vector<double> adaptive_penalty_memory;
  std::vector<double> adaptive_penalty_retention;
  std::list<ThreeCSpike> spikes;
  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;
  double ns_fractional_improvement_final,
      ns_fractional_improvement_state_final, ns_ridge_beta;
  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;
  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;
  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;
 
  void initialize_inverse_operator();
  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;
  double compute_ns_peak_threshold();
  void rescale_spike(ThreeCSpike &spk);
  void update_residual_matrix(const ThreeCSpike &spk);
  void update_lag_weights(const int col, const double candidate_amplitude);
};
} // namespace mspass::algorithms::deconvolution
#endif