Fix decon serialization (#589)

* Fix mistakes in boost serialization

* Fix a number of problems mostly with copy constructors

* Add bindings for ComplexArray and ShapingWavelet classes

* Enable pickle tests for decon processing classes

* Fix debug residual left in my mistake

* Relax prediction error size constraint to reduce random failures

* split pickle dumps and loads in test script

* Add debug code to see why test is failing on github

* Fix mistake using local path for pf file

* Add additional messages to debug serialization problem on github

* Add print statemens to FFTDeconOperator for github debug

* Test possible patch to serialization code

* Add missing attribute and more debug print statements

* Disable serialization scaffolding and pickle test

* Change noise level to make test failure nearly impossible

cxx/include/mspass/algorithms/deconvolution/CNRDeconEngine.h

@@ -153,29 +153,39 @@ class CNRDeconEngine : public FFTDeconOperator
   template<class Archive>
   void serialize(Archive& ar, const unsigned int version)
   {
-      /* See boost documentation for why these are needed.
-      tapers are referenced through a polymorphic shared_ptr.
-      It seems all subclasses to be used need this incantation*/
-      ar.register_type(static_cast<LinearTaper *>(NULL));
-      ar.register_type(static_cast<CosineTaper *>(NULL));
-      ar.register_type(static_cast<VectorTaper *>(NULL));
+      //std::cout <<"Entered serialize function"<<std::endl;
       ar & boost::serialization::base_object<FFTDeconOperator>(*this);
+      //std::cout << "Serializing first group of simple parameters"<<std::endl;
       ar & algorithm;
       ar & damp;
       ar & noise_floor;
       ar & band_snr_floor;
       ar & operator_dt;
       ar & shaping_wavelet_number_poles;
+      //std::cout << "Serializing shapingwavelet"<<std::endl;
       ar & shapingwavelet;
+      ar & winlength;
+      //std::cout<<"Serializing power spectrum engine objects"<<std::endl;
       ar & signal_engine;
       ar & noise_engine;
       ar & snr_regularization_floor;
+      //std::cout << "Serializin winv vector"<<std::endl;
       ar & winv;
       ar & winv_t0_lag;
+      //std::cout<<"Serializing final block of parameters"<<std::endl;
       ar & regularization_bandwidth_fraction;
-      ar & peak_snr;
-      ar & signal_bandwidth_fraction;
-
+      /* These fixed length arrays caused probems - seg faults.
+       * Apparently boost doesn't handle that corectly.  There may 
+       * be a more concise way to do this but this should always work. */
+      //std::cout << "Entering block for 3 component arrays"<<std::endl;
+      for(auto k=0;k<3;++k)
+      {
+        //std::cout << "k="<<k<<std::endl;
+        ar & peak_snr[k];
+        //std::cout << "bandwidth_fraction"<<std::endl;
+        ar & signal_bandwidth_fraction[k];
+      }
+      //std::cout << "Exiting serialize function"<<std::endl;
   }
 };
 }  // End namespace enscapsulation

cxx/include/mspass/algorithms/deconvolution/FFTDeconOperator.h

@@ -64,18 +64,25 @@ class FFTDeconOperator
     template<class Archive>
     void save(Archive& ar, const unsigned int version) const
     {
+      //std::cout << "Entered FFTDecon serialization save function"<<std::endl;
         ar & nfft;
         ar & sample_shift;
         ar & winv;
+        //std::cout << "Exiting save function" << std::endl;
     }
     template<class Archive>
     void load(Archive &ar, const unsigned int version)
     {
+        //std::cout << "Entered FFTDecon serializaton load function" << std::endl;
       ar & this->nfft;
       ar & this->sample_shift;
+        //std::cout << "Loading winv vector" << std::endl;
       ar & winv;
+        //std::cout << "Creating wavetable " << std::endl;
       this->wavetable = gsl_fft_complex_wavetable_alloc (this->nfft);
+        //std::cout << "Creating workspace" << std::endl;
       this->workspace = gsl_fft_complex_workspace_alloc (this->nfft);
+        //std::cout << "Exiting load" << std::endl;
     }
     BOOST_SERIALIZATION_SPLIT_MEMBER()
 };

cxx/include/mspass/algorithms/deconvolution/MultiTaperSpecDivDecon.h

@@ -13,11 +13,11 @@
 #include "mspass/algorithms/deconvolution/ShapingWavelet.h"
 #include "mspass/seismic/CoreTimeSeries.h"
 namespace mspass::algorithms::deconvolution{
-class MultiTaperSpecDivDecon: public ScalarDecon, public FFTDeconOperator
+class MultiTaperSpecDivDecon: public FFTDeconOperator, public ScalarDecon
 {
 public:
     /*! Default constructor.   Do not use - only for declarations */
-    MultiTaperSpecDivDecon() : ScalarDecon(),FFTDeconOperator(){};
+    MultiTaperSpecDivDecon() : FFTDeconOperator(),ScalarDecon(){};
     MultiTaperSpecDivDecon(const mspass::utility::Metadata &md,const std::vector<double> &noise,
                           const std::vector<double> &wavelet,const std::vector<double> &data);
     MultiTaperSpecDivDecon(const mspass::utility::Metadata &md);
@@ -105,17 +105,16 @@ class MultiTaperSpecDivDecon: public ScalarDecon, public FFTDeconOperator
         return nw;
     };
 private:
-    /*! Private method called by constructors to load parameters.   */
-    int read_metadata(const mspass::utility::Metadata &md,bool refresh);
-    /* Returns a tapered data in container of ComplexArray objects*/
-    std::vector<ComplexArray> taper_data(const std::vector<double>& signal);
     std::vector<double> noise;
     double nw,damp;
     int nseq;  // number of tapers
     unsigned int taperlen;
     mspass::utility::dmatrix tapers;
-    /* inverse wavelet in the frequency domain */
-    std::vector<ComplexArray> winv;
+    /* With this algorithm we have to keep frequeny domain 
+     * representations of the inverse for each taper.  The 
+     * xcor version of this alorithm doesn't need that because 
+     * the averaging is different. */
+    std::vector<ComplexArray> winv_taper;
     /* We also cache the actual output fft because the cost is
      * small compared to need to recompute it when requested.
      * This is a feature added for the GID method that adds
@@ -125,20 +124,25 @@ class MultiTaperSpecDivDecon: public ScalarDecon, public FFTDeconOperator
     are averaged to produce final result, but we keep them for the
     option of bootstrap errors. */
     std::vector<ComplexArray> rfestimates;
+    /* Private methods */
+    /* Returns a tapered data in container of ComplexArray objects*/
+    std::vector<ComplexArray> taper_data(const std::vector<double>& signal);
+    /*! Private method called by constructors to load parameters.   */
+    int read_metadata(const mspass::utility::Metadata &md,bool refresh);
     int apply();
     friend boost::serialization::access;
     template<class Archive>
     void serialize(Archive &ar, const unsigned int version)
     {
-        ar & boost::serialization::base_object<FFTDeconOperator>(*this);
         ar & boost::serialization::base_object<ScalarDecon>(*this);
+        ar & boost::serialization::base_object<FFTDeconOperator>(*this);
         ar & noise;
         ar & nw;
         ar & damp;
         ar & nseq;
         ar & taperlen;
         ar & tapers;
-        ar & winv;
+        ar & winv_taper;
         ar & ao_fft;
         ar & rfestimates;
     }

cxx/include/mspass/algorithms/deconvolution/MultiTaperXcorDecon.h

@@ -98,10 +98,7 @@ class MultiTaperXcorDecon: public FFTDeconOperator, public ScalarDecon
         return nw;
     };
 private:
-    /*! Private method called by constructors to load parameters.   */
-    int read_metadata(const mspass::utility::Metadata &md,bool refresh);
-    /* Returns a tapered data in container of ComplexArray objects*/
-    std::vector<ComplexArray> taper_data(const std::vector<double>& signal);
+    /* noise data vector */
     std::vector<double> noise;
     double nw,damp;
     int nseq;  // number of tapers
@@ -112,20 +109,23 @@ class MultiTaperXcorDecon: public FFTDeconOperator, public ScalarDecon
      * This is a feature added for the GID method that adds
      * an inefficiency for straight application */
     ComplexArray ao_fft;
+    /*! Private method called by constructors to load parameters.   */
+    int read_metadata(const mspass::utility::Metadata &md,bool refresh);
+    /* Returns a tapered data in container of ComplexArray objects*/
+    std::vector<ComplexArray> taper_data(const std::vector<double>& signal);
     int apply();
     friend boost::serialization::access;
     template<class Archive>
     void serialize(Archive &ar, const unsigned int version)
     {
-        ar & boost::serialization::base_object<FFTDeconOperator>(*this);
         ar & boost::serialization::base_object<ScalarDecon>(*this);
+        ar & boost::serialization::base_object<FFTDeconOperator>(*this);
         ar & noise;
         ar & nw;
         ar & damp;
         ar & nseq;
         ar & taperlen;
         ar & tapers;
-        ar & winv;
         ar & ao_fft;
     }
 };

cxx/python/algorithms/deconvolution_py.cc

@@ -7,6 +7,8 @@
 
 #include <boost/archive/text_oarchive.hpp>
 
+#include <mspass/algorithms/deconvolution/ComplexArray.h>
+#include <mspass/algorithms/deconvolution/ShapingWavelet.h>
 #include <mspass/algorithms/deconvolution/WaterLevelDecon.h>
 #include <mspass/algorithms/deconvolution/LeastSquareDecon.h>
 #include <mspass/algorithms/deconvolution/MultiTaperXcorDecon.h>
@@ -93,6 +95,35 @@ PYBIND11_MODULE(deconvolution, m) {
 
   /* Need this to support returns of std::vector in children of ScalarDecon*/
   py::bind_vector<std::vector<double>>(m, "DoubleVector");
+  /* All the frequency domain operators use this class internally. 
+   * Useful still to have bindings to the underlying class. */
+  py::class_<ComplexArray>(m,"ComplexArray","Complex valued Fortran style array implementation used in MsPASS Decon opertors")
+    //.def(py::init<std::vector<Complex64&>())
+    .def(py::init<int,double *>())
+    .def(py::init<const ComplexArray&>())
+    .def("conj",&ComplexArray::conj,"Convert array elements to complex conjugates")
+    .def("abs",&ComplexArray::abs,"Return DoubleVector of complex magnitudes")
+    .def("rms",&ComplexArray::rms,"Return rms of array content")
+    .def("norm2",&ComplexArray::norm2,"Return L2 norm of array content")
+    .def("phase",&ComplexArray::phase,"Return DoubleVector of phase of components")
+    .def("size",&ComplexArray::size,"Return number of components in the array")
+    ;
+  /* All frequency domain methods uses this class internally as well. 
+   * It actually contains a ComplexArray.   Useful to have these bindings 
+   * for testing and inspection of the result of a get_shaping_wavelet method.
+   * */
+  py::class_<ShapingWavelet>(m,"ShapingWavelet","Shaping wavelet object used in MsPASS decon frequency domain decon operators")
+    .def(py::init<const Metadata&,int>())
+    .def(py::init<int,double,int,double,double,int>())
+    .def(py::init<double,double,int>())
+    .def("impulse_response",&ShapingWavelet::impulse_response,"Return the impulse response of the wavelet in a CoreTimeSeries container")
+    .def("df",&ShapingWavelet::freq_bin_size,"Return frequency bin size (Hz)")
+    .def("dt",&ShapingWavelet::sample_interval,"Return sample interval of wavelet in the time domain")
+    .def("type",&ShapingWavelet::type,"Return the string description of the type of signal this wavelet defines")
+    .def("size",&ShapingWavelet::size,"Size of the complex array defining the wavelet internally")
+    ;
+
+
 
   /* this is a set of deconvolution related classes*/
   py::class_<ScalarDecon,PyScalarDecon>(m,"ScalarDecon","Base class for scalar TimeSeries data")
@@ -158,6 +189,8 @@ PYBIND11_MODULE(deconvolution, m) {
         boost::archive::text_iarchive artm(sstm);
         LeastSquareDecon lsd;
         artm >> lsd;
+        ShapingWavelet sw=lsd.get_shaping_wavelet();
+        ComplexArray w(*sw.wavelet());
         return lsd;
       }
     ))

cxx/src/lib/algorithms/deconvolution/CNRDeconEngine.cc

@@ -132,7 +132,8 @@ CNRDeconEngine::CNRDeconEngine(const AntelopePf& pf)
 }
 /* Standard copy constructor */
 CNRDeconEngine::CNRDeconEngine(const CNRDeconEngine& parent)
-  : shapingwavelet(parent.shapingwavelet),
+  : FFTDeconOperator(parent),
+     shapingwavelet(parent.shapingwavelet),
       signal_engine(parent.signal_engine),
         noise_engine(parent.noise_engine),
           winv(parent.winv)
@@ -142,6 +143,7 @@ CNRDeconEngine::CNRDeconEngine(const CNRDeconEngine& parent)
   this->noise_floor = parent.noise_floor;
   this->band_snr_floor = parent.band_snr_floor;
   this->operator_dt = parent.operator_dt;
+  this->winlength = parent.winlength;
   this->shaping_wavelet_number_poles = parent.shaping_wavelet_number_poles;
   this->snr_regularization_floor = parent.snr_regularization_floor;
   this->regularization_bandwidth_fraction = parent.regularization_bandwidth_fraction;

cxx/src/lib/algorithms/deconvolution/ComplexArray.cc

@@ -98,7 +98,6 @@ ComplexArray::ComplexArray(vector<double> mag,vector<double> phase)
     }
     else
     {
-        cout<<"Length of magnitude vector and phase vector doesn't match"<<endl;
         throw MsPASSError("ComplexArray::ComplexArray(vector<double> mag,vector<double> phase): Length of magnitude vector and phase vector do not match",
               ErrorSeverity::Invalid);
     }

cxx/src/lib/algorithms/deconvolution/LeastSquareDecon.cc

@@ -10,7 +10,7 @@ using namespace mspass::utility;
 using mspass::algorithms::amplitudes::normalize;
 
 LeastSquareDecon::LeastSquareDecon(const LeastSquareDecon &parent)
-    : FFTDeconOperator(parent)
+    : FFTDeconOperator(parent),ScalarDecon(parent)
 {
     damp=parent.damp;
 }

cxx/src/lib/algorithms/deconvolution/MultiTaperSpecDivDecon.cc

@@ -18,7 +18,7 @@ using namespace mspass::utility;
 using mspass::algorithms::amplitudes::normalize;
 
 MultiTaperSpecDivDecon::MultiTaperSpecDivDecon(const Metadata &md)
-    : ScalarDecon(md), FFTDeconOperator(md)
+    : FFTDeconOperator(md), ScalarDecon(md)
 {
 
     try {
@@ -34,17 +34,18 @@ MultiTaperSpecDivDecon::MultiTaperSpecDivDecon(const Metadata &md)
 }
 
 MultiTaperSpecDivDecon::MultiTaperSpecDivDecon(const MultiTaperSpecDivDecon &parent)
-    : ScalarDecon(parent), FFTDeconOperator(parent), tapers(parent.tapers)
+    : FFTDeconOperator(parent), 
+        ScalarDecon(parent), 
+          tapers(parent.tapers),
+            noise(parent.noise),
+              ao_fft(parent.ao_fft),
+                rfestimates(parent.rfestimates),
+                  winv_taper(parent.winv_taper)
 {
-    /* wavelet and data vectors are copied in ScalarDecon copy constructor.
-    This method needs a noise vector so we have explicitly copy it here. */
-    noise=parent.noise;
-    /* ditto for shaping wavelet vector */
-    shapingwavelet=parent.shapingwavelet;
-    /* multitaper parameters to copy */
     nw=parent.nw;
-    taperlen=parent.taperlen;
+    nseq=parent.nseq;
     damp=parent.damp;
+    taperlen=parent.taperlen;
 }
 int MultiTaperSpecDivDecon::read_metadata(const Metadata &md,bool refresh)
 {
@@ -307,7 +308,7 @@ void MultiTaperSpecDivDecon::process()
     }
     /* Probably should save these in private area for this estimator*/
     //vector<ComplexArray> rfestimates;
-    /* Must clear this and winv containers or they accumulate */
+    /* Must clear rfestimate and winv_taper containers or they accumulate */
     rfestimates.clear();
     for(i=0;i<nseq;++i)
     {
@@ -319,8 +320,7 @@ void MultiTaperSpecDivDecon::process()
     For consistency with related methods we'll store the frequency domain
     values, BUT these now become essentially a matrix - actually stored
     as a vector of vectors */
-    //ComplexArray winv;
-    winv.clear();
+    winv_taper.clear();
     ao_fft.clear();
     double *d0=new double[nfft];
     for(int k=0;k<nfft;++k) d0[k]=0.0;
@@ -332,7 +332,7 @@ void MultiTaperSpecDivDecon::process()
     {
       ComplexArray work(delta0);
       work=work/denominator[i];
-      winv.push_back(work);
+      winv_taper.push_back(work);
     }
     for(i=0; i<nseq; ++i)
     {
@@ -422,7 +422,7 @@ CoreTimeSeries MultiTaperSpecDivDecon::inverse_wavelet(const double t0parent)
      CoreTimeSeries result(this->nfft);
      for(int i=0;i<nseq;++i)
      {
-       CoreTimeSeries work(this->FFTDeconOperator::FourierInverse(this->winv[i],
+       CoreTimeSeries work(this->FFTDeconOperator::FourierInverse(this->winv_taper[i],
                   *shapingwavelet.wavelet(),dt,t0parent));
        if(i==0)
 	   result=work;
@@ -453,7 +453,7 @@ std::vector<CoreTimeSeries> MultiTaperSpecDivDecon::all_inverse_wavelets
         for(int i=0;i<nseq;++i)
         {
             CoreTimeSeries work(this->FFTDeconOperator::FourierInverse
-               (this->winv[i],*shapingwavelet.wavelet(),dt,t0parent));
+               (this->winv_taper[i],*shapingwavelet.wavelet(),dt,t0parent));
             all.push_back(work);
         }
         return all;

cxx/src/lib/algorithms/deconvolution/MultiTaperXcorDecon.cc

@@ -32,15 +32,15 @@ MultiTaperXcorDecon::MultiTaperXcorDecon(const Metadata &md)
 }
 
 MultiTaperXcorDecon::MultiTaperXcorDecon(const MultiTaperXcorDecon &parent)
-    : FFTDeconOperator(parent), tapers(parent.tapers)
+    : FFTDeconOperator(parent), 
+        ScalarDecon(parent), 
+          tapers(parent.tapers),
+            noise(parent.noise),
+              ao_fft(parent.ao_fft)
 {
-    /* wavelet and data vectors are copied in ScalarDecon copy constructor.
-    This method needs a noise vector so we have explicitly copy it here. */
-    noise=parent.noise;
-    /* ditto for shaping wavelet vector */
-    shapingwavelet=parent.shapingwavelet;
     /* multitaper parameters to copy */
     nw=parent.nw;
+    nseq=parent.nseq;
     taperlen=parent.taperlen;
     damp=parent.damp;
 }

cxx/src/lib/algorithms/deconvolution/ScalarDecon.cc

@@ -23,7 +23,10 @@ ScalarDecon::ScalarDecon(const vector<double>& d, const vector<double>& w)
     result.reserve(data.size());
 }
 ScalarDecon::ScalarDecon(const ScalarDecon& parent)
-    : data(parent.data),wavelet(parent.wavelet),result(parent.result)
+    : data(parent.data),
+      wavelet(parent.wavelet),
+       result(parent.result),
+         shapingwavelet(parent.shapingwavelet)
 {
 }
 ScalarDecon& ScalarDecon::operator=(const ScalarDecon& parent)

cxx/src/lib/algorithms/deconvolution/WaterLevelDecon.cc

@@ -9,7 +9,7 @@ using namespace mspass::utility;
 using mspass::algorithms::amplitudes::normalize;
 
 WaterLevelDecon::WaterLevelDecon(const WaterLevelDecon &parent)
-    : FFTDeconOperator(parent)
+    : FFTDeconOperator(parent),ScalarDecon(parent)
 {
     wlv=parent.wlv;
 }