ComplexArray.h

#ifndef __COMPLEX_ARRAY_H__
#define __COMPLEX_ARRAY_H__
 
#include <boost/archive/text_iarchive.hpp>
#include <boost/archive/text_oarchive.hpp>
#include <boost/serialization/shared_ptr.hpp>
#include <complex>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>
#include <iostream>
#include <vector>
 
#define REAL(z, i) ((z)[2 * (i)])
#define IMAG(z, i) ((z)[2 * (i) + 1])
namespace mspass::algorithms::deconvolution {
typedef std::complex<double> Complex64;
typedef std::complex<float> Complex32;
 
typedef struct FortranComplex32 {
  float real;
  float imag;
} FortranComplex32;
typedef struct FortranComplex64 {
  double real;
  double imag;
} FortranComplex64;
/* needed for boost serialization */
template <class Archive>
void serialize(Archive &ar, FortranComplex64 &z, const unsigned int version) {
  ar & z.real;
  ar & z.imag;
}
 
/* \brief Interfacing object to ease conversion between FORTRAN and C++ complex.
 
   */
class ComplexArray {
public:
  ComplexArray();
  ComplexArray(std::vector<Complex64> &d);
  ComplexArray(std::vector<Complex32> &d);
  ComplexArray(int nsamp, FortranComplex32 *d);
  ComplexArray(int nsamp, FortranComplex64 *d);
  ComplexArray(int nsamp, float *d);
  ComplexArray(int nsamp, double *d);
  ComplexArray(int nsamp);
  template <class T> ComplexArray(int nsamp, std::vector<T> d);
  template <class T> ComplexArray(int nsamp, T d);
 
  ComplexArray(std::vector<double> mag, std::vector<double> phase);
 
  /* These will need to be implemented.  Likely cannot
     depend on the compiler to generate them correctly */
  ComplexArray(const ComplexArray &parent);
  ComplexArray &operator=(const ComplexArray &parent);
  ~ComplexArray();
 
  /* These are kind of the inverse of the constructor.
  Independent of what the internal representation is they
  will return useful interface representations. */
  // template<class T> T *FortranData();
  /* This is same for what I think fortran calls
     double complex */
  //        double *FortranData();
  /* C representation.  This can be templated easily.
  See below.  The syntax is weird and should probably
  be wrapped with a typedef */
  // template<class T> std::vector<std::complex<T> > CPPData();
 
  /* Operators are the most important elements of this
     thing to make life easier. */
  Complex64 operator[](int sample);
  double *ptr();
  double *ptr(int sample);
  ComplexArray &operator+=(const ComplexArray &other) noexcept(false);
  ComplexArray &operator-=(const ComplexArray &other) noexcept(false);
  /* This actually is like .* in matlab - sample by sample multiply not
  a dot product */
  ComplexArray &operator*=(const ComplexArray &other) noexcept(false);
  /* This is like *= but complex divide element by element */
  ComplexArray &operator/=(const ComplexArray &other) noexcept(false);
  const ComplexArray operator+(const ComplexArray &other) const noexcept(false);
  // template<class T> ComplexArray operator +(const vector<T> &other);
  // template<class T> ComplexArray operator +(const T &other);
  const ComplexArray operator-(const ComplexArray &other) const noexcept(false);
  // template<class T> ComplexArray operator -(const vector<T> &other);
  // template<class T> ComplexArray operator -(const T &other);
  const ComplexArray operator*(const ComplexArray &other) const noexcept(false);
  const ComplexArray operator/(const ComplexArray &other) const noexcept(false);
  // template<class T> ComplexArray operator *(const vector<T> &other);
  // template<class T> friend ComplexArray operator *(const vector<T> &lhs,const
  // ComplexArray &rhs);
  // template<class T> ComplexArray operator *(const T &other);
  // template<class T> friend ComplexArray operator *(const T &lhs,const
  // ComplexArray &rhs);
  void conj();
  /* Return stl vector of amplitude spectrum.  */
  std::vector<double> abs() const;
  /* Return rms value.  */
  double rms() const;
  /* Return 2-norm value.  */
  double norm2() const;
  /* Return stl vector of phase */
  std::vector<double> phase() const;
  /* Return size of the array*/
  int size() const;
 
private:
  /* Here is an implementation detail.   There are three ways
     I can think to do this.  First, we could internally store
     data as fortran array of 32 bit floats.   That is probably
     the best because we can use BLAS routines (if you haven't
     heard of this - likely - I need to educate you.)  to do
     most of the numerics fast. Second, we could use stl
     vector container of std::complex.  The third is excessively
     messy but technically feasible - I would not recommend it.
     That is, one could store pointers to either representation
     and internally convert back and forth.  Ugly and dangerous
     I think.
 
     I suggest we store a FORTRAN 32 bit form since that is
     what standard numeric libraries (e.g. most fft routines)
     use.  */
  /*I decided to use 64 bit, since the GSL's fft routine is using that.*/
  /* Changed from raw pointer to shared_ptr by glp - Dec 2024 */
  // FortranComplex64 *data;
  std::shared_ptr<FortranComplex64[]> data;
  int nsamp;
  friend boost::serialization::access;
  template <class Archive>
  void save(Archive &ar, const unsigned int version) const {
    std::vector<FortranComplex64> dv;
    dv.reserve(this->nsamp);
    for (auto i = 0; i < nsamp; ++i)
      dv.push_back(this->data[i]);
    ar & nsamp;
    ar & dv;
  }
  template <class Archive> void load(Archive &ar, const unsigned int version) {
    ar &this->nsamp;
    std::vector<FortranComplex64> dv;
    dv.reserve(this->nsamp);
    ar & dv;
    this->data =
        std::shared_ptr<FortranComplex64[]>(new FortranComplex64[this->nsamp]);
    for (auto i = 0; i < this->nsamp; ++i)
      this->data[i] = dv[i];
  }
  BOOST_SERIALIZATION_SPLIT_MEMBER()
};
/* This would normally be in the .h file and since I don't think
   you've used templates worth showing you how it would work. */
/*
template <class T> std::vector<std::complex<T> > ComplexArray::CPPData()
{
    std::vector<std::complex<T> > result;
    result.reserve(nsamp);
    std::size_t  i;
    for(i=0; i<nsamp; ++i)
    {
        std::complex<T> z(data[i].real, data[i].imag);
        result.push_back(z);
    }
    return result;
}
*/
 
/*
template<class T> T* ComplexArray::FortranData()
{
    T* result=new T[nsamp];
    for(std::size_t i=0; i<nsamp; i++)
        result[i]=data[i];
    return result;
}
*/
 
template <class T> ComplexArray::ComplexArray(int n, std::vector<T> d) {
  nsamp = n;
  if (nsamp > d.size()) {
    data = std::shared_ptr<FortranComplex64[]>(new FortranComplex64[nsamp]);
    for (std::size_t i = 0; i < d.size(); i++) {
      this->data[i].real = d[i];
      this->data[i].imag = 0.0;
    }
    for (std::size_t i = d.size(); i < nsamp; i++) {
      this->data[i].real = 0.0;
      this->data[i].imag = 0.0;
    }
  } else {
    data = std::shared_ptr<FortranComplex64[]>(new FortranComplex64[nsamp]);
    for (std::size_t i = 0; i < nsamp; i++) {
      this->data[i].real = d[i];
      this->data[i].imag = 0.0;
    }
  }
}
template <class T> ComplexArray::ComplexArray(int n, T d) {
  nsamp = n;
  data = std::shared_ptr<FortranComplex64[]>(new FortranComplex64[nsamp]);
  for (std::size_t i = 0; i < nsamp; i++) {
    this->data[i].real = d;
    this->data[i].imag = 0.0;
  }
}
/*
template<class T> ComplexArray ComplexArray::operator +(const vector<T> &other)
{
    ComplexArray result(*this);
    int n;
    if(nsamp>other.size())
        n=other.size();
    else
        n=nsamp;
    for(int i=0; i<n; i++)
    {
        result.data[i].real=data[i].real+other[i];
    }
    return result;
}
template<class T> ComplexArray ComplexArray::operator +(const T &other)
{
    ComplexArray result(*this);
    for(int i=0; i<nsamp; i++)
    {
        result.data[i].real=data[i].real+other;
    }
    return result;
}
template<class T> ComplexArray ComplexArray::operator -(const vector<T> &other)
{
    ComplexArray result(*this);
    int n;
    if(nsamp>other.size())
        n=other.size();
    else
        n=nsamp;
    for(int i=0; i<n; i++)
    {
        result.data[i].real=data[i].real-other[i];
    }
    return result;
}
template<class T> ComplexArray ComplexArray::operator -(const T &other)
{
    ComplexArray result(*this);
    for(int i=0; i<nsamp; i++)
    {
        result.data[i].real=data[i].real-other;
    }
    return result;
}
template<class T> ComplexArray ComplexArray::operator *(const vector<T> &other)
{
    ComplexArray result(*this);
    int n;
    if(nsamp>other.size())
        n=other.size();
    else
        n=nsamp;
    for(int i=0; i<n; i++)
    {
        result.data[i].real=data[i].real*other[i];
        result.data[i].imag=data[i].imag*other[i];
    }
    return result;
}
template<class T> ComplexArray operator *(const vector<T>& lhs,const
ComplexArray& rhs)
{
    return rhs*lhs;
}
template<class T> ComplexArray ComplexArray::operator *(const T &other)
{
    ComplexArray result(*this);
    for(int i=0; i<nsamp; i++)
    {
        result.data[i].real=data[i].real*other;
        result.data[i].imag=data[i].imag*other;
    }
    return result;
}
template<class T> ComplexArray operator *(const T& lhs,const ComplexArray& rhs)
{
    return rhs*lhs;
}
*/
} // namespace mspass::algorithms::deconvolution
#endif