psc_solver.hpp

/* psc_solver.hpp */
/* $Id$ */
#ifndef __psc_solver_hpp
#define __psc_solver_hpp

#include <utility>
#include <vector>
#include <set>
#include <string>
#include <lapackpp/laslv.h>
#include <lapackpp/laexcp.h>

#ifdef GE_BUILD
#include <wave_field.hpp>
#else
#include <s2proc/wave_field.hpp>
#endif

namespace psc
{
  class Solver
  {
  public:
    
    enum RowType
      {
    FIRST_ORDER_ROW_TYPE = 0, 
    SECOND_ORDER_ROW_TYPE     
      };
    
    struct TimeShift
    {
      int id;          
      Trace::TimeT dt; 
    };
    
    struct RowElem
    {
      int id1, id2, id3, id4; 
      Trace::TimeT dt;        
      static bool stride_asc_id1_asc_comp(const RowElem &re1, const RowElem &re2) 
      {
    int stride1 = abs(re1.id2 - re1.id1);
    int stride2 = abs(re2.id2 - re2.id1);
    
    if (stride1 < stride2)
      return true;
    else if (stride1 > stride2)
      return false;
    else 
      return (re1.id1 < re2.id1);
      }
    };
    
    typedef std::vector<RowElem> RowElemVec;

    typedef std::vector<TimeShift> TimeShifts; 
    
    struct Params
    {
      Params():
    row_type(FIRST_ORDER_ROW_TYPE),
    max_time_shift(-1),
    max_used_stride(0),
    remove_trend_flag(true),
    trend_polynomial_degree(1), 
    remove_low_freq_flag(false),
    low_freq_polynomial_degree(1),
    ids_range(std::pair<int,int>(0,0))
      {/* nothing to do */}
      
      RowType      row_type;                   
      Trace::TimeT max_time_shift;             
      unsigned int max_used_stride;            
      bool         remove_trend_flag;          
      size_t       trend_polynomial_degree;    
      bool         remove_low_freq_flag;       
      size_t       low_freq_polynomial_degree; 
      std::pair<int,int> ids_range;            
    };
    
    
    Solver(const Params &params);
    
    virtual ~Solver();
    
    
    bool operator()(const WaveField &row, 
            TimeShifts &time_shifts) const;
    
    bool operator()(RowElemVec &row_elem_vec, 
            TimeShifts &time_shifts) const;
    
    void setParameters(const Params &params);
    
    const Params& getParameters() const;
    
    
  protected:
        
    Params m_params;
    
    
    void update_eq_data(const WaveField &row, 
            RowElemVec &row_elem_vec) const; 
    
    bool determine_statics(const RowElemVec &row_elem_vec, 
               TimeShifts &time_shifts) const;
    
    std::set<int> get_row_ids(const RowElemVec &row_elem_vec) const;
    
    bool get_eq_row(const RowElem &row_elem, const std::set<int> &ids_set, 
            std::vector<Trace::AmplT> &row_vec, Trace::TimeT &rp_elem) const;
    
    bool remove_low_freq(const size_t polynomial_degree, RowElemVec &row_elem_vec) const;
    
    bool remove_low_freq(const size_t polynomial_degree, 
             RowElemVec::iterator elem_it1,
             RowElemVec::iterator elem_it2) const;
    
    bool remove_trend(const size_t trend_polynomial_degree, TimeShifts& time_shifts) const;
    
    
    void test_la_to_files(const std::string &matr_fname, const std::string &vec_fname, 
              const LaGenMatDouble &Am, const LaVectorDouble &bv) const;    
  };
}; // namespace psc

#endif /* psc_solver.hpp */