ConvolveWithConst.cxx

Go to the documentation of this file.

#include "CgemDigitizerSvc/ConvolveWithConst.h"
 
//#include <cstdint>
#include <complex>
#include <vector>
#include <cstring>
#include "TVirtualFFT.h"
#include <climits>
#include <iostream>
 
using std::cout;
using std::cerr;
using std::endl;
typedef unsigned int uint32_tt;
//static_assert( UINT_MAX == 4294967295U," a unsigned 32-bit int is needed for Bit Twidding Hack at ConvolveWithConst::ConvolveWithConst() ");
#if not UINT_MAX == 4294967295U // i.e. uint is uint32_t 
    Compile Time Error A Unsigned 32-bit int is needed  at ConvolveWithConst::ConvolveWithConst;
#endif
 
static const char options1[4][10]={"R2C ES K","R2C M K","R2C P K","R2C EX K"};
static const char options2[4][10]={"C2R ES K","C2R M K","C2R P K","C2R EX K"};
 
 
struct ConvolveWithConst::convParams{
    int length;
    ConvolveWithConst::TransformOptimOpt optOptim;
    //bool isR2C;// not ifft.
    bool operator< (const convParams b) const{return (long)(this->length)*4+(long)(this->optOptim) <
                                (long)(b.length)*4+(long)(b.optOptim);// long should be longer than int here?
    }
};
struct ConvolveWithConst::fft2{
    TVirtualFFT * fft;
    TVirtualFFT * ifft;
    unsigned int refcount;// ++ when CWC init, -- when CWC final, delete when ==0U and CWC final
    fft2(){fft=nullptr;ifft=nullptr;refcount=0U;}
 
};
std::map<ConvolveWithConst::convParams,ConvolveWithConst::fft2> ConvolveWithConst::m_map;
//const std::string ConvolveWithConst::options1[1]="R2C ES K";
//const std::string ConvolveWithConst::options1[2]="R2C M K" ;
//const std::string ConvolveWithConst::options1[3]="R2C P K" ;
//const std::string ConvolveWithConst::options1[4]="R2C EX K";
//const std::string ConvolveWithConst::options2[1]="C2R ES K";
//const std::string ConvolveWithConst::options2[2]="C2R M K" ;
//const std::string ConvolveWithConst::options2[3]="C2R P K" ;
//const std::string ConvolveWithConst::options2[4]="C2R EX K";
 
ConvolveWithConst::ConvolveWithConst(const double* ConstArray,
            const int ConstArrayLength,const int Blength, 
            TransformOptimOpt opt,TransformLengthOpt optL){
            m_initialized=false;
            m_fft=nullptr;
            m_ifft=nullptr;
            this->init( ConstArray,
            ConstArrayLength,Blength, 
            opt, optL);
            }
 
void ConvolveWithConst::init(const double* ConstArray,
            const int ConstArrayLength,const int Blength, 
            TransformOptimOpt opt,TransformLengthOpt optL){
    if (m_initialized) {
        std::cout << "CgemDigitizerSvc::ConvolveWithConst::init(): trying to initialize a instance that claimed itself initialized. may cause problem!!!!" << std::endl;
        return;}// dont initialize it too many times. 
    m_optOptim=opt;
    m_optLength=optL;
    m_BLength=Blength;
    switch(m_optLength){
        uint32_tt len;
        case ConvolveWithConst::opt_Nearst2Pow:
            len=ConstArrayLength+Blength-1;
            // Bit Twidding Hack. basically gets a greatest greater 2 power
            len--; 
            len |= len >> 1;
            len |= len >> 2;
            len |= len >> 4;
            len |= len >> 8;
            len |= len >> 16;
            m_length=len+1;
            break;
        case ConvolveWithConst::opt_AsShortAsPsb:
            m_length=ConstArrayLength+Blength-1;
            break;
    }
    std::vector<double> A_tmp;
 
    convParams pars={m_length,opt};
    if (0==m_map[pars].fft){
        m_map[pars].fft =TVirtualFFT::FFT(1, &m_length,options1[opt]);
        m_map[pars].ifft=TVirtualFFT::FFT(1, &m_length,options2[opt]);
        
        //std::cout<<"newed: fft at "<<m_map[pars].fft<<std::endl;
    }
    m_fft  = m_map[pars].fft;
    m_ifft = m_map[pars].ifft;
    m_ref  = m_map.find(pars);// should be same... as []
    m_map[pars].refcount++;
    
    //m_fft=TVirtualFFT::FFT(1, &m_length,"R2C EX K");
    //m_ifft=TVirtualFFT::FFT(1, &m_length,"C2R EX K");
 
 
    //m_fft=TVirtualFFT::FFT(1, &m_length,options1[m_optOptim]);
    //m_ifft=TVirtualFFT::FFT(1, &m_length,options2[m_optOptim]);
    zeros.clear();
    zeros.resize(m_length*2);
    A_fft.resize(m_length*2);
 
    A_tmp.clear();
    
    A_tmp.resize(m_length); 
    
    
    double *p_A_fft=&(A_fft.front()); 
    double *p_A_tmp=&(A_tmp.front()); 
 
    std::memcpy(p_A_tmp,ConstArray,ConstArrayLength*sizeof(double));
    m_fft->SetPoints(p_A_tmp);
    m_fft->Transform();
    m_fft->GetPoints(p_A_fft);
 
    //m_length=-1;
    m_initialized=true;
 
}
ConvolveWithConst::ConvolveWithConst(){
    m_initialized=false;
    m_fft=nullptr;
    m_ifft=nullptr;
}
 
 
ConvolveWithConst::~ConvolveWithConst(){
    //std::cout<<"m_fft:"<<m_fft<<"  m_ifft:"<<m_ifft<<std::endl;
    //sometimes due to problem of TVirtualFFT? the TVirtualFFT* is the global fft.
    // consider rewrite. now we may roughly check if it is on stack 
    //if ((size_t)m_fft < 0x700000000000LU)
    //    delete m_fft;
    //if ((size_t)m_ifft < 0x700000000000LU)
    //    delete m_ifft;
    //delete[] A_fft, B_fft,AB_fft;//,out;
    //std::cout<<"ConvolveWithConst::~ConvolveWithConst(), m_fft="<<m_fft<<" ";
    if (not (m_fft==nullptr)){
        m_ref->second.refcount--;
        //std::cout<<"refcount --; now is "<<m_ref->second.refcount<<" ;";
        if (m_ref->second.refcount==0){
            //std::cout<<"thus deleted m_fft="<<m_ref->second.fft<<" ";
            delete m_ref->second.fft;
            delete m_ref->second.ifft;
            m_map.erase(m_ref);
        }
    }
    //std::cout<<std::endl;
 
}
/**
 * @brief do a convolve of stored const A and B, and put results to output.
 * 
 * @param output address of output, array of double needed since memcpy
 * @param B array of input B
 * @param leftindex output is copied from resuts[leftindex]
 * @param sizeofOut to results[leftindex+sizeofOut], assuming results has indefinite trailing 0
 * @param factor each value is multiplicated with factor, making output[:]*=factor
 */
void ConvolveWithConst::convolve(double * output , 
            const double *  B, 
            int leftIndex,int sizeofOut,double factor)const{
    if (not m_initialized){throw "ConvolveWithConst::convolve: object not initialized";}
    if (leftIndex > m_length) throw "ConvolveWithConst::convolve: leftIndex" 
                "should not be more than NA+NB-1, as there should be only 0";
    if (leftIndex <0 ) throw "ConvolveWithConst::convolve: leftIndex should not be negative.";
    //if (leftIndex + sizeofOut > zeros.size()+m_length) throw "ConvolveWithConst::convolve: leftIndex + sizeofOut too much larger than  NA+NB-1";
    B_tmp.clear();
    B_tmp.resize(m_length);
    B_fft.resize(m_length*2);
    AB_fft.resize(m_length*2);
 
    double *p_A_fft=&(A_fft.front()); 
    double *p_B_fft=&(B_fft.front()); 
    double *p_AB_fft=&(AB_fft.front());
    //double *p_A_tmp=&(A_tmp.front()); 
    double *p_B_tmp=&(B_tmp.front()); 
    const double *p_zeros=&(zeros.front()); 
    
    //std::memcpy(p_A_tmp,A,sizeofA*sizeof(double));// we do clear and memcpy to add "padding" zeros.
    std::memcpy(p_B_tmp,B,m_BLength*sizeof(double));
 
    m_fft->SetPoints(p_B_tmp);
    m_fft->Transform();
    m_fft->GetPoints(p_B_fft);
    // the 2 fft just writes the left part, 
    //leaving the right part for us to mirror and conj.
    //leaving a optimize point.
 
    for (int i = 0; i < (m_length / 2 + 1); i++) {
        ((Complex *)p_AB_fft)[i] = ((Complex *)p_A_fft)[i] * ((Complex *)p_B_fft)[i] / (double)m_length * factor; //simple complex product. 
    }
    for (int i = m_length - 1; i >= (m_length / 2 + 1); i--) {
        AB_fft[2 * i]     = AB_fft[2 * (m_length - i)];
        AB_fft[2 * i + 1] = -AB_fft[2 * (m_length - i) + 1];// mirror and conj
    }
    m_ifft->SetPoints(p_AB_fft);
    m_ifft->Transform();
    m_ifft->GetPoints(p_B_tmp);// since B_tmp already used and no longer needed in this session, we reuse this vector.
    if (leftIndex+sizeofOut <=m_length){
        std::memcpy(output,(p_B_tmp+leftIndex),sizeofOut*sizeof(double));
    }
    else{// add padding 0
        std::memcpy(output,(p_B_tmp+leftIndex),(m_length-leftIndex)*sizeof(double));
        while (leftIndex+sizeofOut -m_length > zeros.size()){
            std::memcpy(output+m_length-leftIndex,(p_zeros),(zeros.size())*sizeof(double));
            output+=zeros.size();
            sizeofOut-=zeros.size();
 
        }
        std::memcpy(output+m_length-leftIndex,(p_zeros),(leftIndex+sizeofOut -m_length)*sizeof(double));
    }
 
}
/**
 * @brief multiply with ffted saved results and add them to output with factor; 
 *          out += in* SavedConst*factor/L
 * @param out 
 * @param in 
 * @param factor 
 */
void ConvolveWithConst::MultiplyAndAdd(ConvolveWithConst::Complex* outHalfPlus1, const ConvolveWithConst::Complex* inHalfPlus1, double factor) const{
    double *p_A_fft=&(A_fft.front()); 
    for (int i = 0; i < (m_length / 2 + 1); i++) {
        outHalfPlus1[i] += inHalfPlus1[i] * ((Complex *)p_A_fft)[i] / (double)m_length * factor; //simple complex product. 
    }
}
//void ConvolveWithConst::destoryAll(){
//    for(std::map<convParams,fft2>::iterator it=m_map.begin();it !=m_map.end();it++){
//        
//            delete it->second.fft;
//            delete it->second.ifft;
//    }
//    m_map.clear();
//}
 
/**
 * @brief Do fft with respect of getLength(); 
 * 
 * @param outHalfPlus1 
 * @param inNormalLength 
 * @param LengthOutAsDouble 
 * @param LengthIn 
 */
void ConvolveWithConst::FFT(ConvolveWithConst::Complex* outHalfPlus1, const double * inNormalLength, int LengthOutAsDouble, int LengthIn) const{
    if (not m_initialized){cerr<<"ConvolveWithConst::FFT: object not initialized"<<endl; throw "ConvolveWithConst::FFT: object not initialized";}
    if (LengthOutAsDouble/2 < m_length / 2 +1) {cerr<<"ConvolveWithConst::FFT:LengthOutAsDouble too small" <<endl;throw "LengthOutAsDouble too small";} // opt_Nearst2Pow tested only.
    B_tmp.clear();
    B_tmp.resize(m_length);
    double *p_B_tmp=&(B_tmp.front()); 
    std::memcpy(p_B_tmp,inNormalLength,LengthIn*sizeof(double));
 
    m_fft->SetPoints(p_B_tmp);
    m_fft->Transform();
    m_fft->GetPoints((double*)outHalfPlus1);
}
/**
 * @brief Do ifft with respect of getLength(); 
 * 
 * @param outHalfPlus1 
 * @param inNormalLength 
 * @param LengthOutAsDouble 
 * @param LengthIn 
 */
void ConvolveWithConst::IFFT(double* outNormalLength, const ConvolveWithConst::Complex * inHalfPlus1, int lengthOut, int lengthInAsDouble) const{
    if (not m_initialized){cerr<<"ConvolveWithConst::IFFT: object not initialized"<<endl; throw "ConvolveWithConst::IFFT: object not initialized";}
    if (lengthInAsDouble/2 < m_length/2+1){ cerr<<"ConvolveWithConst::IFFT:LengthInAsDouble too small"<<endl;throw "LengthInAsDouble too small";}
    AB_fft.resize(m_length*2);
    double *p_AB_fft=&(AB_fft.front());
    B_tmp.clear();
    B_tmp.resize(m_length);
    double *p_B_tmp=&(B_tmp.front()); 
 
    //for (int i = 0; i < (m_length / 2 + 1); i++) {
    //    ((Complex *)p_AB_fft)[i] = ((Complex *)p_A_fft)[i] * ((Complex *)p_B_fft)[i] / (double)m_length * factor; //simple complex product. 
    //}
    std::memcpy(p_AB_fft,inHalfPlus1,lengthInAsDouble*sizeof(double));
    for (int i = m_length - 1; i >= (m_length / 2 + 1); i--) {
        AB_fft[2 * i]     = AB_fft[2 * (m_length - i)];
        AB_fft[2 * i + 1] = -AB_fft[2 * (m_length - i) + 1];// mirror and conj
    }
 
 
 
 
    m_ifft->SetPoints(p_AB_fft);
    m_ifft->Transform();
    m_ifft->GetPoints(p_B_tmp);
    if (lengthOut> m_length)lengthOut = m_length;
    std::memcpy(outNormalLength,p_B_tmp,lengthOut*sizeof(double));
 
}
//void ConvolveWithConst::mirrorAndConjugateInPlace(ConvolveWithConst::Complex * inFullLength,int lengthInAsDouble) const{
//    if (lengthInAsDouble/2 < length) throw std::invalid_argument("LengthInAsDouble too small");
//
//    for (int i = m_length - 1; i >= (m_length / 2 + 1); i--) {
//        inFullLength[2 * i]     = inFullLength[2 * (m_length - i)];
//        inFullLength[2 * i + 1] = -inFullLength[2 * (m_length - i) + 1];// mirror and conj
//    }
//}