QVectorC.cc

Go to the documentation of this file.

/*
* 
* Class QVectorC
*
*/
 
#include "QVectorC.hh"
#include "QMatrixC.hh"
#include "QError.hh"
#include <math.h>
#include <sstream>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_vector_complex_double.h>
#include <gsl/gsl_complex_math.h>
 
#ifdef _HAVE_FFTW3_
#include <fftw3.h>
#endif
 
QObjectImp(Diana::QVectorC);
 
Q_BEGIN_NAMESPACE
 
 
void QVectorC::ArrayCopy(const gsl_complex* orig,const UInt_t origStride, 
                         gsl_complex *dest,const UInt_t destStride,
                         const UInt_t size)
{
//TODO: find faster copy
  if(size > 0 && orig && dest) {
    if(origStride == 1 && destStride == 1)
      memcpy(dest,orig,sizeof(gsl_complex)*size);
    else
      for(UInt_t i =0; i < size; i++) {
        dest[i*destStride] = orig[i*origStride]; 
      }
  }
}
QVectorC::QVectorC()
{
  fDataOwner = true;
  fData=0;
  fSize = 0;
  fMathVec = 0;
  fStride = 1;
  fAllocSize = 0;
  fRe=NULL;
  fIm=NULL;
  fInterleaved=NULL;
}
 
QVectorC::QVectorC(const UInt_t size)
{
  fDataOwner = true;
  fData = QVector::ArrayAlloc(size*2);  
  fAllocSize = size*2;
  fSize = size*2;
  fMathVec = 0;
  fStride = 1;
  fRe=NULL;
  fIm=NULL;
  fInterleaved=NULL;
}
QVectorC::QVectorC(const QVectorC& in)
{
  InitCopy(in,true);
}
QVectorC::QVectorC(const QVector& vec)
{  
  fSize=(vec.Size());
  fAllocSize=fSize;
  fDataOwner = true;
  fData = QVector::ArrayAlloc(fSize);  
  fMathVec = NULL;
  fStride = 1;
  for (UInt_t i=0; i<Size(); i++)  {
    fData[2*i] = vec[i];
    fData[2*i+1]= vec[Size()+i];
  }
  fRe=NULL;
  fIm=NULL;
  fInterleaved=NULL;
}
QVectorC::QVectorC(const QVector& re, const QVector& im)
{    
 
  if(re.Size()==im.Size())
    {
      fSize=re.fSize*2;
      fAllocSize=fSize;
      fData = QVector::ArrayAlloc(fSize);
      fStride=1;  
      for (UInt_t i=0; i<Size(); i++)
        {
          fData[2*i]=re[i];
          fData[2*i+1]=im[i];
        }
      fDataOwner = true;
      fMathVec = NULL;
      fRe=NULL;
      fIm=NULL;
      fInterleaved=NULL;
    }else{
    DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"You are trying to define a QVectorC  with different sizes QVectors") )
  }
}
 
QVectorC::QVectorC(const QVectorC& orig,const bool dataOwner)
{
  InitCopy(orig,dataOwner);
}
QVectorC::QVectorC(gsl_complex *data,const UInt_t size,const UInt_t stride){
  fData=(double *)data;
  fDataOwner = false;
  fStride=stride;
  fSize = 2*size;
  fAllocSize = (UInt_t)(2*size*stride);// should never been used
  fMathVec = NULL;
  fRe=NULL;
  fIm=NULL;
  fInterleaved=NULL;
}
void QVectorC::InitCopy(const QVectorC& orig,const bool dataOwner)
{
  fDataOwner = dataOwner;       
  if(fDataOwner) {    
    fStride=1;
    fData = 0;
    fSize = 0;
    fAllocSize = 0;
    if(orig.fSize!=0 && orig.fData ) {
      fSize = orig.fSize;
      fData = QVector::ArrayAlloc(fSize);
      fAllocSize = fSize;
      // This copies complex numbers, so we need fSize/2
      ArrayCopy((gsl_complex *)orig.fData,orig.fStride,
                (gsl_complex *)fData,1,fSize/2);
    }
  }
  else {    
    fStride=orig.fStride;
    fData = orig.fData;
    fSize = orig.fSize;
    fAllocSize = orig.fAllocSize;
  }
  fMathVec=NULL;
  fRe=NULL;
  fIm=NULL;
  fInterleaved=NULL;
}
QVectorC::~QVectorC()
{
  if(fDataOwner) {
    if(fData) QVector::ArrayFree(fData);
  }
  if(fMathVec) {
    free(fMathVec);
  }
  if(fRe) {
    delete fRe;
  }
  if(fIm) {
    delete fIm;
  }
  if(fInterleaved){
    delete fInterleaved;
  }
  fSize=0;
  fAllocSize=0;
  fData=NULL;
  fMathVec=NULL;
  fRe=NULL;
  fIm=NULL;
  fInterleaved=NULL;
}
 
void QVectorC::Resize(const UInt_t newSize)
{
  if(fDataOwner) {
    if(newSize*2 > fAllocSize) {
      if(fData) QVector::ArrayFree(fData);
      fAllocSize=newSize*2;
      fData = QVector::ArrayAlloc(fAllocSize);
    }
    fSize = newSize*2;
  } 
  else if(newSize != Size() ) {
    DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"Not data owner, are you trying to resize a QVectorC get from a QVectorC(Const)View?") );
  }
  if(fMathVec) free(fMathVec);
  fMathVec = NULL;
  if(fRe) {
    delete fRe;
  }
  if(fIm) {
    delete fIm;
  }
  if(fInterleaved) {
    delete fInterleaved;
  }
  fRe=NULL;
  fIm=NULL;
  fInterleaved=NULL;
}
 
void QVectorC::SetArray(const double* re, const double* im,const UInt_t size) 
{
  Resize(size);
  for(UInt_t i=0;i<size;i++){
    fData[i*2*fStride]=re[i];
    fData[i*2*fStride+1]=im[i];
  }
}
void QVectorC::Initialize(const double re, const double im)
{
  for(UInt_t i=0;i<Size();i++) {
    fData[2*i*fStride]=re;
    fData[2*i*fStride+1]=im;
  }
}
void QVectorC::Initialize(const QComplex val){
  QVectorC::Initialize(val.Re(),val.Im());
}
 
const QVectorC& QVectorC::Rebin(const int rebin,bool sum) {
  int newSize = Size()/rebin;
  double norm = (sum ? 1. : 1./rebin);
  for(int i=0;i<newSize;i++){
    (*this)[i]=norm*(*this)[rebin*i];
    for(int j=1;j<rebin;j++){
      (*this)[i]+=norm*(*this)[rebin*i+j];
    }
  }
  Resize(newSize);
  return *this;
}
 
const QVectorC& QVectorC::Shift(const int nstep)
{
  bool gtz = nstep > 0 ? true : false;
  if(nstep == 0 || Size() == 0) return *this;
  unsigned int unstep = abs(nstep) % Size();
  QVectorC rest(unstep);
  if(gtz) {
    for(unsigned int i = 0; i < unstep; i++) rest[i] = (*this)[Size() - unstep + i];
    for(int i = Size() - unstep - 1; i >= 0 ; i--) (*this)[i + unstep] = (*this)[i];
    for(unsigned int i = 0; i < unstep; i++) (*this)[i] = rest[i];
  } else {      
    for(unsigned int i = 0; i < unstep; i++) rest[i] = (*this)[i];
    for(unsigned int i = unstep; i < Size(); i++) (*this)[i - unstep] = (*this)[i];
    for(unsigned int i = 0; i < unstep; i++) (*this)[i + Size() - unstep] = rest[i];
  }
  return *this;
}
 
const QVectorC& QVectorC::ShiftReal(const double fstep)
{
  int nstep = (int) ceil(fstep);
  Shift(nstep);
  double fracs = nstep - fstep;
  Diana::QComplex zero = (*this)[0];
  for(int i = 0; i < (int)Size() - 1; i++) {
    (*this)[i] = (*this)[i]*(1.-fracs) + fracs*(*this)[i+1]; 
  }
  (*this)[Size()-1] = (*this)[Size()-1]*(1.-fracs) + fracs*zero;
  
  return *this;
}
 
QComplex QVectorC::operator[](const UInt_t i)
{
  if(i < Size() ) {
    return (fData+i*2*fStride);
  }
  std::stringstream stream;
  stream<<"Size: "<<Size()<<" Index: "<<i;
  DianaThrow( QError(QERR_OUT_OF_RANGE,__FILE__,__LINE__,stream.str()) )
  return Q_DOUBLE_DEFAULT;
}
 
const QComplex QVectorC::operator[](UInt_t i) const
{
  if(i < Size() ) {
    return (fData+i*2*fStride);
  }
  std::stringstream stream;
  stream<<"Size: "<<Size()<<" Index: "<<i;
  DianaThrow( QError(QERR_OUT_OF_RANGE,__FILE__,__LINE__,stream.str()) );
  return Q_DOUBLE_DEFAULT;
}
 
const QVectorC& QVectorC::operator=(const QVectorC& orig) 
{
  Resize( orig.Size() );
  ArrayCopy((gsl_complex *)orig.fData,orig.fStride,
            (gsl_complex *)fData,fStride,fSize/2);
  return *this;
}
QMatrixC QVectorC::operator*(const QMatrixC& matrix) const
{
  QMatrixC tmp(Size(),1);
  tmp.SetCol(0,*this);
  tmp *= matrix;
  return tmp;
  }
QVector& QVectorC::Re() 
{
  SetParts();
  return *fRe;
}
QVector& QVectorC::Im() 
{
  SetParts();
  return *fIm;
}
const QVector& QVectorC::Re() const 
{
  SetParts();
  return *fRe;
}
const QVector& QVectorC::Im() const 
{
  SetParts();
  return *fIm;
}
QVector QVectorC::GetPhase() const
{
  QVector res(Size());
  for(UInt_t i = 0; i < Size(); i++) {
    res[i] = operator[](i).GetPhase();
  }
  return res;
}
 
QVector QVectorC::GetMagnitudes() const
{
  QVector res(Size());
  for(UInt_t i=0;i<Size();i++){
    res[i]=sqrt( fData[i*2*fStride]*fData[i*2*fStride] + fData[i*2*fStride+1]*fData[i*2*fStride+1] );
  }
  return res;
}
QVector QVectorC::GetMagnitudesSquare() const
{
  QVector res(Size());
  for(UInt_t i=0;i<Size();i++) {
    res[i]=fData[i*2*fStride]*fData[i*2*fStride] + fData[i*2*fStride+1]*fData[i*2*fStride+1];
  }
  return res;
}
 
QVectorC QVectorC::operator-() const 
{
  QVectorC tmp(*this);
  tmp*=(-1.);
  return tmp;
}
QVectorC QVectorC::Mult(const QVectorC& b) const
{
  //solo con le gsl 1.14
  QVectorC tmp(*this);
  //  tmp.Initialize();
  if(b.Size()==Size()){
    b.SetMathVector();
    tmp.SetMathVector();
    gsl_vector_complex_mul(tmp.fMathVec, b.fMathVec);
  }else{
    DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"You are trying to use QVectorC::Mult with different sizes QVectorCs") );
  }
  return tmp;
}
const QVectorC& QVectorC::MultIn(const QVectorC& b)
{
  if(b.Size()==Size()){
    SetMathVector();
    b.SetMathVector();
    gsl_vector_complex_mul(fMathVec, b.fMathVec);
  }else{
    DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"You are trying to use QVectorC::Mult with different sizes QVectorCs") );
  }
  return *this;
}
 
QVectorC QVectorC::Mult(const QVector& b) const
{
  //solo con le gsl 1.14
  QVectorC tmp(*this);
  //TODO: is there a good gsl method for this?
  if(b.Size()==Size()){
    for (UInt_t i=0; i<Size(); i++){
      tmp[i]*=b[i];
    }
    //SetMathVector();
    //tmp.SetMathVector();
    //gsl_vector_complex_mul(tmp.fMathVec, b.fMathVec);
  }else{
    DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"You are trying to use QVectorC::Mult with different sizes QVectorCs") );
  }
  return tmp;
}
QVectorC QVectorC::Div(const QVectorC&b) const 
{
  QVectorC tmp(*this);
  if(b.fSize==fSize) {
    b.SetMathVector();
    tmp.SetMathVector();
    gsl_vector_complex_div(tmp.fMathVec, b.fMathVec);
  }else{
    DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"You are trying to use QVectorC::Div with of different sizes QVectorCs"));
  }
  return tmp; 
}
 
const QVectorC& QVectorC::DivIn(const QVectorC&b)
{
  if(b.fSize==fSize) {
    SetMathVector();
    b.SetMathVector();
    gsl_vector_complex_div(fMathVec, b.fMathVec);
  }else{
    DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"You are trying to use QVectorC::Div with of different sizes QVectorCs"));
  }
  return *this; 
}
QVectorC QVectorC::Div(const QVector &b) const 
{
  QVectorC tmp(*this);
  //TODO: is there a good gsl method for this?
  if(b.Size()==(UInt_t)Size()) {
   for (UInt_t i=0;i<Size(); i++) {
     tmp[i]/=b[i];
   }
//     SetMathVector();
//     tmp.SetMathVector();
//     gsl_vector_complex_div(tmp.fMathVec, b.fMathVec);
  }else{
    DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"You are trying to use QVectorC::Div with of different sizes QVectorCs"));
  }
  return tmp; 
}
 
const QVectorC& QVectorC::operator/=(const double t) 
{
  return operator*=(1./t);
}
 
const QVectorC& QVectorC::operator*=(const double t) 
{
  SetMathVector();
  gsl_complex scale = gsl_complex_rect(t,0);
  gsl_vector_complex_scale(fMathVec, scale);
  return *this;
}
 
const QVectorC& QVectorC::operator+=(const QVectorC& b) 
{
  if(b.fSize==fSize) {
//     for (UInt_t i=0; i<Size(); i++){
//       fData[2*i*fStride]+=b.fData[2*i*b.fStride];
//       fData[2*i*fStride+1]+=b.fData[2*i*b.fStride+1];
//     }
    b.SetMathVector();
    SetMathVector();
    gsl_vector_complex_add(fMathVec, b.fMathVec);
  }else{
    DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"You are trying to use QVectorC::operator+= with different sizes QVectorCs?"));
  }
  return *this;
}
 
const QVectorC& QVectorC::operator-=(const QVectorC& other) 
{
  operator+=(other*(-1));
  return *this;
}
 
QComplex QVectorC::operator*(const QVectorC& other) const
{
  gsl_complex res;
  other.SetMathVector();
  SetMathVector();
  gsl_blas_zdotu(fMathVec,other.fMathVec,&res);
  return QComplex(res);
}
 
QVectorC QVectorC::operator*(const double t) const
{ 
  QVectorC tmp(*this);
  tmp*=t;
  return tmp;
}
 
QVectorC QVectorC::operator/(const double t) const
{
  QVectorC tmp(*this);
  tmp/=t;
  return tmp;
}
 
QVectorC QVectorC::operator+(const QVectorC& other) const
{
  QVectorC tmp(*this);
  tmp+=other;
  return tmp;
}
 
QVectorC QVectorC::operator-(const QVectorC& other) const
{
  QVectorC tmp(*this);
  tmp-=other;
  return tmp;
}
 
const QVectorC& QVectorC::operator*=(const QComplex& t)
{
  SetMathVector();
  gsl_vector_complex_scale(fMathVec, t.Get_gsl_complex());
//   for(UInt_t i=0;i<Size();i++){
//     operator[](i)*=t;
//   }
  return *this;
}
const QVectorC& QVectorC::operator-=(const QComplex& t)
{
    gsl_complex sub = gsl_complex_mul_real(t.Get_gsl_complex(), -1);
    SetMathVector();
    gsl_vector_complex_add_constant(fMathVec, sub);
//   for(UInt_t i=0;i<Size();i++){
//     operator[](i)-=t;
//   }
   return *this;
}
 
const QVectorC& QVectorC::operator/=(const QComplex& t)
{
  SetMathVector();
  gsl_complex inv = gsl_complex_inverse(t.Get_gsl_complex());
  gsl_vector_complex_scale(fMathVec, inv);
//   for(UInt_t i=0;i<Size();i++)  {
//     operator[](i)/=t;
//   }
  return *this;
}
QVectorC QVectorC::operator*(const QComplex& t) const
{
  QVectorC tmp(*this);
  tmp*=t;
  return tmp;
}
 
QVectorC QVectorC::operator/(const QComplex& t) const
{
  QVectorC tmp(*this);
  tmp/=t;
  return tmp;
}
QMatrixC QVectorC::H() const
{
  QMatrixC tmp(1,Size());
  tmp.SetRow(0,Conj());
  return tmp;
}
QVectorC &QVectorC::Conjugate()
{
  //TODO: Is there a gls function for this?
  for(UInt_t i=0;i<Size();i++)
    fData[i*2*fStride+1]*=-1;
  return *this;
}
QVectorC QVectorC::Conj() const
{
  QVectorC tmp(*this);
  tmp.Conjugate();
  return tmp;
}
 
void QVectorC::Append(const QComplex& val)
{
  if(fDataOwner){
    double* temp = QVector::ArrayAlloc((Size()+1)*2);
    ArrayCopy((gsl_complex *)fData,fStride,(gsl_complex *)temp,1,Size());
    if(fData) QVector::ArrayFree(fData);
    fSize=(Size()+1)*2;
    fAllocSize = fSize;;
    fStride=1;//check
    temp[fSize-2]=val.Re();
    temp[fSize-1]=val.Im();
    fData=temp;
    fStride=1;
    if(fMathVec) free(fMathVec);
    fMathVec = NULL;
  }
  else {
    DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                       "Not data owner, you are trying to resize a "
                       "QVectorC using Append"));
  }
}
void QVectorC::Append(const double re,const double im)
{
  QComplex val(re,im);
  Append(val);
  return;
}
void QVectorC::Append(const QVectorC& vec)
{
  for(UInt_t i=0;i<vec.Size();i++) {
    QComplex val=vec[i];
    Append(val);
  }
  return;
}
void QVectorC::SetRe(const QVector& re)
{  
  if(re.Size()== Size() ){
    for (UInt_t i=0;i<Size(); i++) {
      fData[2*i*fStride]=re[i];
    }
    return;
  }
  std::stringstream msg;
  msg<<"You are trying to use QVectorC::SetRe with a wrong size QVector: ("
     <<re.Size()<<") instead of ("<<Size()<<")"<<std::endl;
  DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,msg.str()));
}
void QVectorC::SetIm(const QVector& im)
{  
  if(im.Size()==Size() ) {
    for (UInt_t i=0; i <Size(); i++)  {
      fData[2*i*fStride+1]=im[i];
    }
    return;
  }
  std::stringstream msg;
  msg<<"You are trying to use QVectorC::SetIm with a wrong size QVector: ("
          <<im.Size()<<") instead of ("<<Size()<<")"<<std::endl;
 DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,msg.str()))
}
double QVectorC::Norma() const
{
  SetMathVector();
  return gsl_blas_dznrm2(fMathVec);
}
QVector QVectorC::SingleVector() const
{
  QVector res(fSize);
  for (UInt_t i=0;i<Size();i++) {
    res[i]=fData[2*i*fStride];
    res[i+Size()]=fData[2*i*fStride+1];
  }
  return res;
}
QVector &QVectorC::InterleavedVector() const 
{
  if(fStride!=1)
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      "Can not interleave QVectorC if Stride !=1!"));
  if(fInterleaved==NULL) 
    fInterleaved=new QVector(fData,fSize,1);
  return *fInterleaved;
}
QVectorC QVectorC::GetSubVector(UInt_t N,UInt_t start,UInt_t stride){
  if(start+N > Size())
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      "QVectorC: Sub-vector would overrun parent vector!"));
  stride = (stride < 1 ? 1 : stride);
  QVectorC ret((gsl_complex *)fData+start,N,fStride*stride);
  return ret;
}
const QVectorC QVectorC::GetSubVector(UInt_t N,UInt_t start,UInt_t stride)
  const {
  if(start+N > Size())
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      "QVectorC: Sub-vector would overrun parent vector!"));
  stride = (stride < 1 ? 1 : stride);
  QVectorC ret((gsl_complex *)fData+start,N,fStride*stride);
  return ret;
}
 
QComplex QVectorC::Sum(const UInt_t nelem,const UInt_t first) const
{
  QComplex sum(0,0);
  UInt_t last = first + nelem;
  for(UInt_t i = first; i < last; i++) sum += this->operator[](i);
  return sum;
}
 
 
void QVectorC::SetMathVector() const
{
  if(!fMathVec) {
    fMathVec = (gsl_vector_complex*) malloc(sizeof(gsl_vector_complex));      
    fMathVec->size = Size();
    fMathVec->stride = fStride;
    fMathVec->data = fData;
    fMathVec->block = 0;
    fMathVec->owner = 0;
  }
}
void QVectorC::SetParts() const
{
  if(!fRe)  fRe = new QVector(fData,Size(),fStride*2);  
  if(!fIm)  fIm = new QVector(fData+1,Size(),fStride*2);
}
void QVectorC::Print(const char *opts) const {
 
  std::cout << "fSize " << Size() << std::endl;
  std::cout << "fStride " << fStride << std::endl;
  std::cout << "fData " << fData << std::endl;
  std::cout << "fData   ";
  for(UInt_t i=0;i<Size();i++)
    std::cout << (*this)[i].Re() << " " << (*this)[i].Im() << " ";
  std::cout << std::endl;
 
}
Q_END_NAMESPACE
 
/********************************************
  NON MEMBER FUNCTIONS
 *********************************************/
 
std::ostream& operator<<(std::ostream&s,const Diana::QVectorC &vec)
{
  for (UInt_t i = 0; i <(UInt_t) vec.Size(); i++) s<<vec[i]<<" ";
  s<<std::endl;
  return s;
}
 
Diana::QVectorC operator*(double t, const Diana::QVectorC &v)
{
  Diana::QVectorC tmp(v);
  tmp*=t;
  return tmp;
}
const Diana::QVector Re(const Diana::QVectorC &Z) { return Z.Re(); }
 
const Diana::QVector Im(const Diana::QVectorC &Z) { return Z.Im(); }