CGEM BOSS 6.6.5.g
BESIII Offline Software System
Loading...
Searching...
No Matches
Classes | Public Types | Public Member Functions | List of all members
ConvolveWithConst Class Reference

basically, convolve with fft. convolves 1-d. if size of 2 inputs are fixed, it may yield optimal speed. currently not thread safe. More...

#include <ConvolveWithConst.h>

Classes

struct  convParams
 
struct  fft2
 

Public Types

enum  TransformOptimOpt { opt_ES , opt_M , opt_P , opt_EX }
 
enum  TransformLengthOpt { opt_AsShortAsPsb , opt_Nearst2Pow }
 
typedef std::complex< double > Complex
 

Public Member Functions

 ConvolveWithConst ()
 
void init (const double *ConstArray, const int ConstArrayLength, const int Blength, TransformOptimOpt opt=opt_EX, TransformLengthOpt optL=opt_AsShortAsPsb)
 
 ConvolveWithConst (const double *ConstArray, const int ConstArrayLength, const int Blength, TransformOptimOpt opt=opt_EX, TransformLengthOpt optL=opt_AsShortAsPsb)
 
 ~ConvolveWithConst ()
 
void convolve (double *output, const double *B, const int leftIndex, const int sizeofOut, double factor=1) const
 do a convolve of stored const A and B, and put results to output.
 
int getLength () const
 
void MultiplyAndAdd (Complex *outHalfPlus1, const Complex *inHalfPlus1, double factor=1) const
 multiply with ffted saved results and add them to output with factor; out += in* SavedConst*factor/L
 
void FFT (Complex *outHalfPlus1, const double *inNormalLength, int LengthOutAsDouble, int LengthIn) const
 Do fft with respect of getLength();.
 
void IFFT (double *outNormalLength, const Complex *inHalfPlus1, int lengthOut, int lengthInAsDouble) const
 Do ifft with respect of getLength();.
 

Detailed Description

basically, convolve with fft. convolves 1-d. if size of 2 inputs are fixed, it may yield optimal speed. currently not thread safe.

Definition at line 18 of file ConvolveWithConst.h.

Member Typedef Documentation

◆ Complex

typedef std::complex<double> ConvolveWithConst::Complex

Definition at line 22 of file ConvolveWithConst.h.

Member Enumeration Documentation

◆ TransformLengthOpt

enum ConvolveWithConst::TransformLengthOpt
Enumerator
opt_AsShortAsPsb 
opt_Nearst2Pow 

Definition at line 21 of file ConvolveWithConst.h.

21{opt_AsShortAsPsb, opt_Nearst2Pow};

◆ TransformOptimOpt

enum ConvolveWithConst::TransformOptimOpt
Enumerator
opt_ES 
opt_M 
opt_P 
opt_EX 

Definition at line 20 of file ConvolveWithConst.h.

20{opt_ES,opt_M,opt_P,opt_EX};

Constructor & Destructor Documentation

◆ ConvolveWithConst() [1/2]

ConvolveWithConst::ConvolveWithConst ( )

Definition at line 127 of file ConvolveWithConst.cxx.

127 {
128 m_initialized=false;
129 m_fft=nullptr;
130 m_ifft=nullptr;
131}

◆ ConvolveWithConst() [2/2]

ConvolveWithConst::ConvolveWithConst ( const double *  ConstArray,
const int  ConstArrayLength,
const int  Blength,
TransformOptimOpt  opt = opt_EX,
TransformLengthOpt  optL = opt_AsShortAsPsb 
)

Definition at line 49 of file ConvolveWithConst.cxx.

51 {
52 m_initialized=false;
53 m_fft=nullptr;
54 m_ifft=nullptr;
55 this->init( ConstArray,
56 ConstArrayLength,Blength,
57 opt, optL);
58 }
ConvolveWithConst::init
void init(const double *ConstArray, const int ConstArrayLength, const int Blength, TransformOptimOpt opt=opt_EX, TransformLengthOpt optL=opt_AsShortAsPsb)
Definition: ConvolveWithConst.cxx:60

◆ ~ConvolveWithConst()

ConvolveWithConst::~ConvolveWithConst ( )

Definition at line 134 of file ConvolveWithConst.cxx.

134 {
135 //std::cout<<"m_fft:"<<m_fft<<" m_ifft:"<<m_ifft<<std::endl;
136 //sometimes due to problem of TVirtualFFT? the TVirtualFFT* is the global fft.
137 // consider rewrite. now we may roughly check if it is on stack
138 //if ((size_t)m_fft < 0x700000000000LU)
139 // delete m_fft;
140 //if ((size_t)m_ifft < 0x700000000000LU)
141 // delete m_ifft;
142 //delete[] A_fft, B_fft,AB_fft;//,out;
143 //std::cout<<"ConvolveWithConst::~ConvolveWithConst(), m_fft="<<m_fft<<" ";
144 if (not (m_fft==nullptr)){
145 m_ref->second.refcount--;
146 //std::cout<<"refcount --; now is "<<m_ref->second.refcount<<" ;";
147 if (m_ref->second.refcount==0){
148 //std::cout<<"thus deleted m_fft="<<m_ref->second.fft<<" ";
149 delete m_ref->second.fft;
150 delete m_ref->second.ifft;
151 m_map.erase(m_ref);
152 }
153 }
154 //std::cout<<std::endl;
155
156}

Member Function Documentation

◆ convolve()

void ConvolveWithConst::convolve ( double *  output,
const double *  B,
const int  leftIndex,
const int  sizeofOut,
double  factor = 1 
) const

do a convolve of stored const A and B, and put results to output.

Parameters
outputaddress of output, array of double needed since memcpy
Barray of input B
leftindexoutput is copied from resuts[leftindex]
sizeofOutto results[leftindex+sizeofOut], assuming results has indefinite trailing 0
factoreach value is multiplicated with factor, making output[:]*=factor

Definition at line 166 of file ConvolveWithConst.cxx.

168 {
169 if (not m_initialized){throw "ConvolveWithConst::convolve: object not initialized";}
170 if (leftIndex > m_length) throw "ConvolveWithConst::convolve: leftIndex"
171 "should not be more than NA+NB-1, as there should be only 0";
172 if (leftIndex <0 ) throw "ConvolveWithConst::convolve: leftIndex should not be negative.";
173 //if (leftIndex + sizeofOut > zeros.size()+m_length) throw "ConvolveWithConst::convolve: leftIndex + sizeofOut too much larger than NA+NB-1";
174 B_tmp.clear();
175 B_tmp.resize(m_length);
176 B_fft.resize(m_length*2);
177 AB_fft.resize(m_length*2);
178
179 double *p_A_fft=&(A_fft.front());
180 double *p_B_fft=&(B_fft.front());
181 double *p_AB_fft=&(AB_fft.front());
182 //double *p_A_tmp=&(A_tmp.front());
183 double *p_B_tmp=&(B_tmp.front());
184 const double *p_zeros=&(zeros.front());
185
186 //std::memcpy(p_A_tmp,A,sizeofA*sizeof(double));// we do clear and memcpy to add "padding" zeros.
187 std::memcpy(p_B_tmp,B,m_BLength*sizeof(double));
188
189 m_fft->SetPoints(p_B_tmp);
190 m_fft->Transform();
191 m_fft->GetPoints(p_B_fft);
192 // the 2 fft just writes the left part,
193 //leaving the right part for us to mirror and conj.
194 //leaving a optimize point.
195
196 for (int i = 0; i < (m_length / 2 + 1); i++) {
197 ((Complex *)p_AB_fft)[i] = ((Complex *)p_A_fft)[i] * ((Complex *)p_B_fft)[i] / (double)m_length * factor; //simple complex product.
198 }
199 for (int i = m_length - 1; i >= (m_length / 2 + 1); i--) {
200 AB_fft[2 * i] = AB_fft[2 * (m_length - i)];
201 AB_fft[2 * i + 1] = -AB_fft[2 * (m_length - i) + 1];// mirror and conj
202 }
203 m_ifft->SetPoints(p_AB_fft);
204 m_ifft->Transform();
205 m_ifft->GetPoints(p_B_tmp);// since B_tmp already used and no longer needed in this session, we reuse this vector.
206 if (leftIndex+sizeofOut <=m_length){
207 std::memcpy(output,(p_B_tmp+leftIndex),sizeofOut*sizeof(double));
208 }
209 else{// add padding 0
210 std::memcpy(output,(p_B_tmp+leftIndex),(m_length-leftIndex)*sizeof(double));
211 while (leftIndex+sizeofOut -m_length > zeros.size()){
212 std::memcpy(output+m_length-leftIndex,(p_zeros),(zeros.size())*sizeof(double));
213 output+=zeros.size();
214 sizeofOut-=zeros.size();
215
216 }
217 std::memcpy(output+m_length-leftIndex,(p_zeros),(leftIndex+sizeofOut -m_length)*sizeof(double));
218 }
219
220}
output
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per saves r n generator level $ !Flag for chat level in output
Definition: FoamA.h:89
ConvolveWithConst::Complex
std::complex< double > Complex
Definition: ConvolveWithConst.h:22

Referenced by CEF::Convolution_Ebranch_fft(), CEF2::Convolution_Ebranch_fft(), CTF::Convolution_Tbranch_fft(), and CTF2::Convolution_Tbranch_fft().

◆ FFT()

void ConvolveWithConst::FFT ( ConvolveWithConst::Complex *  outHalfPlus1,
const double *  inNormalLength,
int  LengthOutAsDouble,
int  LengthIn 
) const

Do fft with respect of getLength();.

Parameters
outHalfPlus1
inNormalLength
LengthOutAsDouble
LengthIn

Definition at line 251 of file ConvolveWithConst.cxx.

251 {
252 if (not m_initialized){cerr<<"ConvolveWithConst::FFT: object not initialized"<<endl; throw "ConvolveWithConst::FFT: object not initialized";}
253 if (LengthOutAsDouble/2 < m_length / 2 +1) {cerr<<"ConvolveWithConst::FFT:LengthOutAsDouble too small" <<endl;throw "LengthOutAsDouble too small";} // opt_Nearst2Pow tested only.
254 B_tmp.clear();
255 B_tmp.resize(m_length);
256 double *p_B_tmp=&(B_tmp.front());
257 std::memcpy(p_B_tmp,inNormalLength,LengthIn*sizeof(double));
258
259 m_fft->SetPoints(p_B_tmp);
260 m_fft->Transform();
261 m_fft->GetPoints((double*)outHalfPlus1);
262}

◆ getLength()

int ConvolveWithConst::getLength ( ) const
inline

Definition at line 55 of file ConvolveWithConst.h.

55{return m_length;};

◆ IFFT()

void ConvolveWithConst::IFFT ( double *  outNormalLength,
const Complex *  inHalfPlus1,
int  lengthOut,
int  lengthInAsDouble 
) const

Do ifft with respect of getLength();.

Parameters
outHalfPlus1
inNormalLength
LengthOutAsDouble
LengthIn

Definition at line 271 of file ConvolveWithConst.cxx.

271 {
272 if (not m_initialized){cerr<<"ConvolveWithConst::IFFT: object not initialized"<<endl; throw "ConvolveWithConst::IFFT: object not initialized";}
273 if (lengthInAsDouble/2 < m_length/2+1){ cerr<<"ConvolveWithConst::IFFT:LengthInAsDouble too small"<<endl;throw "LengthInAsDouble too small";}
274 AB_fft.resize(m_length*2);
275 double *p_AB_fft=&(AB_fft.front());
276 B_tmp.clear();
277 B_tmp.resize(m_length);
278 double *p_B_tmp=&(B_tmp.front());
279
280 //for (int i = 0; i < (m_length / 2 + 1); i++) {
281 // ((Complex *)p_AB_fft)[i] = ((Complex *)p_A_fft)[i] * ((Complex *)p_B_fft)[i] / (double)m_length * factor; //simple complex product.
282 //}
283 std::memcpy(p_AB_fft,inHalfPlus1,lengthInAsDouble*sizeof(double));
284 for (int i = m_length - 1; i >= (m_length / 2 + 1); i--) {
285 AB_fft[2 * i] = AB_fft[2 * (m_length - i)];
286 AB_fft[2 * i + 1] = -AB_fft[2 * (m_length - i) + 1];// mirror and conj
287 }
288
289
290
291
292 m_ifft->SetPoints(p_AB_fft);
293 m_ifft->Transform();
294 m_ifft->GetPoints(p_B_tmp);
295 if (lengthOut> m_length)lengthOut = m_length;
296 std::memcpy(outNormalLength,p_B_tmp,lengthOut*sizeof(double));
297
298}

◆ init()

void ConvolveWithConst::init ( const double *  ConstArray,
const int  ConstArrayLength,
const int  Blength,
TransformOptimOpt  opt = opt_EX,
TransformLengthOpt  optL = opt_AsShortAsPsb 
)

Definition at line 60 of file ConvolveWithConst.cxx.

62 {
63 if (m_initialized) {
64 std::cout << "CgemDigitizerSvc::ConvolveWithConst::init(): trying to initialize a instance that claimed itself initialized. may cause problem!!!!" << std::endl;
65 return;}// dont initialize it too many times.
66 m_optOptim=opt;
67 m_optLength=optL;
68 m_BLength=Blength;
69 switch(m_optLength){
70 uint32_tt len;
71 case ConvolveWithConst::opt_Nearst2Pow:
72 len=ConstArrayLength+Blength-1;
73 // Bit Twidding Hack. basically gets a greatest greater 2 power
74 len--;
75 len |= len >> 1;
76 len |= len >> 2;
77 len |= len >> 4;
78 len |= len >> 8;
79 len |= len >> 16;
80 m_length=len+1;
81 break;
82 case ConvolveWithConst::opt_AsShortAsPsb:
83 m_length=ConstArrayLength+Blength-1;
84 break;
85 }
86 std::vector<double> A_tmp;
87
88 convParams pars={m_length,opt};
89 if (0==m_map[pars].fft){
90 m_map[pars].fft =TVirtualFFT::FFT(1, &m_length,options1[opt]);
91 m_map[pars].ifft=TVirtualFFT::FFT(1, &m_length,options2[opt]);
92
93 //std::cout<<"newed: fft at "<<m_map[pars].fft<<std::endl;
94 }
95 m_fft = m_map[pars].fft;
96 m_ifft = m_map[pars].ifft;
97 m_ref = m_map.find(pars);// should be same... as []
98 m_map[pars].refcount++;
99
100 //m_fft=TVirtualFFT::FFT(1, &m_length,"R2C EX K");
101 //m_ifft=TVirtualFFT::FFT(1, &m_length,"C2R EX K");
102
103
104 //m_fft=TVirtualFFT::FFT(1, &m_length,options1[m_optOptim]);
105 //m_ifft=TVirtualFFT::FFT(1, &m_length,options2[m_optOptim]);
106 zeros.clear();
107 zeros.resize(m_length*2);
108 A_fft.resize(m_length*2);
109
110 A_tmp.clear();
111
112 A_tmp.resize(m_length);
113
114
115 double *p_A_fft=&(A_fft.front());
116 double *p_A_tmp=&(A_tmp.front());
117
118 std::memcpy(p_A_tmp,ConstArray,ConstArrayLength*sizeof(double));
119 m_fft->SetPoints(p_A_tmp);
120 m_fft->Transform();
121 m_fft->GetPoints(p_A_fft);
122
123 //m_length=-1;
124 m_initialized=true;
125
126}
uint32_tt
unsigned int uint32_tt
Definition: ConvolveWithConst.cxx:14

Referenced by ConvolveWithConst(), InductionGar::init(), and InductionGar2::init().

◆ MultiplyAndAdd()

void ConvolveWithConst::MultiplyAndAdd ( ConvolveWithConst::Complex *  outHalfPlus1,
const Complex *  inHalfPlus1,
double  factor = 1 
) const

multiply with ffted saved results and add them to output with factor; out += in* SavedConst*factor/L

Parameters
out
in
factor

Definition at line 228 of file ConvolveWithConst.cxx.

228 {
229 double *p_A_fft=&(A_fft.front());
230 for (int i = 0; i < (m_length / 2 + 1); i++) {
231 outHalfPlus1[i] += inHalfPlus1[i] * ((Complex *)p_A_fft)[i] / (double)m_length * factor; //simple complex product.
232 }
233}
The documentation for this class was generated from the following files: