37 this->
add(in,ch,fL,fH);
56 std::cout<<
"I am =\n";
92 for(
size_t i=1;
i<I;
i++) {
97 this->
chList.push_back(ww);
98 this->
chName.push_back(ch);
99 this->
chMask.push_back(mask);
117 std::cout<<
"regression::add() ERROR: add target first\n";
121 std::cout<<
"regression::add() ERROR: second witness channel can not be added to LPEF\n";
126 std::cout<<
"regression::add() warning: empty witness channel.\n";
135 if(n || abs(m)!=
int(pow(2.,l)+0.1)) {
136 std::cout<<
"regression::add() ERROR: incompatible target and witness channels\n";
162 std::cout<<
"regression::add() ERROR: incompatible target and witness channels\n";
166 ws.
setlow(fL>=fH ? 0. : fL);
172 for(
size_t i=0;
i<I;
i++) {
180 this->
chList.push_back(ws);
181 this->
chName.push_back(ch);
182 this->
chMask.push_back(mask);
201 std::cout<<
"regression::add() ERROR: second witness channel can not be added to LPEF\n";
205 std::cout<<
"regression::add() WARNING: can't construct witness channel.\n";
220 for(
size_t i=0;
i<I;
i++) {
221 if(!(pMn->
data[
i])) {
225 if(!(pMm->
data[
i])) {
229 if(!(pMn->
data[
i]))
continue;
230 for(
size_t j=0;
j<I;
j++) {
231 if(!(pMm->
data[
j]))
continue;
240 for(
size_t i=0;
i<I;
i++) {
247 wsm *= wsn; ts = wsm; ts -= ts.
mean();
249 this->
chList.push_back(wsn);
250 this->
chName.push_back(ch);
251 this->
chMask.push_back(mask);
270 size_t I = this->
chList[0].maxLayer();
276 this->
FILTER.clear(); std::vector<Wiener>().swap(
FILTER);
277 this->
matrix.clear(); std::vector<TMatrixDSym>().swap(
matrix);
278 this->
vCROSS.clear(); std::vector< wavearray<double> >().swap(
vCROSS);
279 this->
vEIGEN.clear(); std::vector< wavearray<double> >().swap(
vEIGEN);
283 if(!pT->
data[i])
continue;
287 for(j=0; j<
chList.size(); j++) {
289 if(!pW->
data[i])
continue;
291 WF.
layer.push_back(i);
292 WF.
norm.push_back(1.);
298 else pT->
data[
i] = 0;
311 int* pM = this->
chMask[
n].data;
314 for(
size_t i=0;
i<I;
i++) {
329 int* pM = this->
chMask[
n].data;
332 for(
size_t i=0;
i<I;
i++) {
349 size_t I = this->
vEIGEN.size();
350 if(n<(
int)I && n>=0)
return vEIGEN[
n];
351 for(
size_t i=0;
i<I;
i++) E.
append(this->vEIGEN[
i]);
372 if(nT>=0 && nT<N) {i=
nT; I=nT+1;}
375 for(
int n=i ;
n<I;
n++) {
378 if(nW>0 && nW<M) {j=nW; J=nW+1;}
381 for(
int m=j;
m<J;
m++) {
386 for(
int k=0;
k<
K;
k++) {
394 OUT.
append((re*re-RE*RE)/(re*re+RE*RE));
417 int I = this->
FILTER.size();
422 int sIZe,sTEp,k00n,k90n,k00m,k90m;
424 double FLTR = !strcmp(this->
chName[0],this->
chName[1]) ? 0. : 1.;
427 double fraction, rms;
441 for(i=0; i<(
int)this->
chList.size(); i++) {
444 if(!I) {std::cout<<
"regression::nothing to clean.\n";
return;}
446 this->
matrix.clear(); std::vector<TMatrixDSym>().swap(
matrix);
447 this->
vCROSS.clear(); std::vector< wavearray<double> >().swap(
vCROSS);
448 this->
vEIGEN.clear(); std::vector< wavearray<double> >().swap(
vEIGEN);
452 N = pW->
layer.size();
456 pTFm = &(this->
chList[0]);
473 for(j=0; j<(
int)ww.
size(); j++) {
476 rms = sqrt(qq.
mean(fm));
488 for(k=-K; k<=
K; k++) {
493 k00m = s00m.
start()-k*sTEp;
494 k90m = s90m.
start()-k*sTEp;
497 k00n = s00n.
start()+k*sTEp;
498 k90n = s90n.
start()+k*sTEp;
502 for(j=K; j<sIZe-
K; j++) {
513 if(k==0) {VC.
data[kk] *= FLTR; VC.
data[kk+K4/2] *= FLTR;}
540 for(k=-K2; k<=K2; k++) {
545 k00m = s00m.
start()-k*sTEp;
546 k90m = s90m.
start()-k*sTEp;
549 k00n = s00n.
start()+k*sTEp;
550 k90n = s90n.
start()+k*sTEp;
555 for(j=K2; j<sIZe-K2; j++) {
570 for(ii=-K; ii<=
K; ii++) {
571 for(jj=-K; jj<=
K; jj++) {
572 MS[nn+ii][mm+jj] = (ii==0 || jj==0) ? acf[ii-jj]*FLTR : acf[ii-jj];
573 MS[mm+jj][nn+ii] = (ii==0 || jj==0) ? acf[ii-jj]*FLTR : acf[ii-jj];
574 MS[nn+ii][mm+jj+K4/2] = (ii==0 || jj==0) ? ccf[ii-jj]*FLTR : ccf[ii-jj];
575 MS[mm+jj+K4/2][nn+ii] = (ii==0 || jj==0) ? ccf[ii-jj]*FLTR : ccf[ii-jj];
576 MS[nn+ii+K4/2][mm+jj] = (ii==0 || jj==0) ? -ccf[ii-jj]*FLTR :-ccf[ii-jj];
577 MS[mm+jj][nn+ii+K4/2] = (ii==0 || jj==0) ? -ccf[ii-jj]*FLTR :-ccf[ii-jj];
578 MS[nn+ii+K4/2][mm+jj+K4/2] = (ii==0 || jj==0) ? acf[ii-jj]*FLTR : acf[ii-jj];
579 MS[mm+jj+K4/2][nn+ii+K4/2] = (ii==0 || jj==0) ? acf[ii-jj]*FLTR : acf[ii-jj];
587 this->
matrix.push_back(MS);
588 this->
vCROSS.push_back(VC);
602 int nA = this->
FILTER.size();
603 int nR = this->
matrix.size();
604 int nC = this->
vCROSS.size();
607 if(!nA) {std::cout<<
"regression::nothing to clean.\n";
return;}
608 if (nA!=nR || nA!=nC || nA==0) {
609 std::cout <<
"Error, filter not initialized\n";
exit(1);
616 this->
vEIGEN.clear(); std::vector< wavearray<double> >().swap(
vEIGEN);
618 for (i=0; i<nA; i++) {
619 int J =
FILTER[
i].layer.size()-1;
620 ne = nE<=0 ? K4*J : nE-1;
621 if(ne>=K4*J) ne = K4*J-1;
627 TMatrixDSymEigen QP(R);
628 TVectorD eigenval(QP.GetEigenValues());
629 TMatrixD eigenvec(QP.GetEigenVectors());
633 double TH = th<0. ? -th*eigenval[0] : th+1.e-12;
635 for(j=0; j<K4*J; j++) {
636 lambda.
data[
j] = eigenval[
j];
637 if(eigenval[j]>=TH) nlast++;
639 if(eigenval[j+1]>eigenval[j])
640 std::cout<<eigenval[
j]<<
" "<<eigenval[j+1]<<endl;
642 if(nlast<1) nlast = 1;
644 this->
vEIGEN.push_back(lambda);
645 if(nlast>ne) nlast = ne;
655 last = eigenval[nlast]>0. ? 1./eigenval[nlast] : 0.;
658 last = 1./eigenval[0];
661 for(j=0;j<=nlast;j++) lambda[j] = eigenval[j]>0. ? 1./eigenval[
j] : 0.;
662 for(j=nlast+1;j<K4*J;j++) lambda[j] = last;
668 for(j=0;j<K4*J;j++) {
670 for(k=0;k<K4*J;k++) temp += eigenvec[k][j]*pC->
data[k];
674 for(j=0;j<K4*J;j++) {
676 for(k=0;k<K4*J;k++) temp += eigenvec[j][k]*vv.
data[k];
681 for (j=0; j<J; j++) {
682 for(k=0;k<=2*
K;k++) {
684 FILTER[
i].filter90[j+1][
k] = aa[j*K4+k+K4/2];
706 int nA = this->
FILTER.size();
707 int nR = this->
matrix.size();
708 int nC = this->
vCROSS.size();
710 if(!nA) {std::cout<<
"regression::nothing to clean.\n";
return;}
711 if (nA!=nR || nA!=nC || nA==0) {
712 std::cout <<
"Error, filter is not initialized\n";
exit(1);
720 for(n=0; n<
L; n++) this->
vrank.push_back(tmp);
735 std::vector< wavearray<double> > wno;
736 std::vector< wavearray<double> > WNO;
749 if(KK<K) KK =
K; KK++;
752 double trms = wno[0].rms(
S);
753 double TRMS = WNO[0].rms(
S);
758 int* pM = this->chMask[j].data;
760 if(c==
'A' || c==
'M') {
761 qq=wno[0]; qq-=wno[
m];
762 QQ=WNO[0]; QQ-=WNO[
m];
765 sum = trms*trms-sum*sum;
766 SUM = TRMS*TRMS-SUM*
SUM;
769 sum = pow(wno[m].rms(
S),2);
770 SUM = pow(WNO[m].rms(
S),2);
778 if(sum+SUM < threshold*threshold) {
779 if(c==
'm' || c==
'M') pM[pF->
layer[
m]] = 0;
782 ww += wno[
m]; WW += WNO[
m];
785 if(c==
'A' || c==
'M') this->
vrank[0].data[
n] = trms*trms+TRMS*TRMS;
798 if(c==
'n') this->
rnoise = wnoise;
799 else if(c==
'N') this->
rnoise = WNOISE;
822 int nA = this->
FILTER.size();
823 int nR = this->
matrix.size();
824 int nC = this->
vCROSS.size();
827 std::cout<<
"regression::nothing to clean.\n";
830 if (nA!=nR || nA!=nC) {
831 std::cout <<
"Error, filter not initialized\n";
exit(1);
839 std::vector< wavearray<double> > vrms;
843 for (
int i=0;
i<tsize;
i++) {
845 if(fL<fH && (freq>fH || freq<fL))
continue;
853 for(m=N-1; m>=N-abs(nbins); m--) {
854 rms.
data[
n] += vrms[
n].data[
m];
855 if(vrms[n].data[m]>0.) nA++;
857 if(nA) rms.
data[
n] /= nA;
879 std::vector<double> *
ff, *FF;
897 pW = &(this->chList[j]);
906 qq *= 1./pF->
norm[
m];
907 QQ *= 1./pF->
norm[
m];
909 for(i=0; i<(
int)qq.
size(); i++) {
911 if(i>=K && i<
int(qq.
size())-K) {
912 for(k=-K; k<=
K; k++) {
913 val += (*ff)[k+
K]*qq.
data[i+
k] - (*FF)[k+
K]*QQ.
data[i+
k];
914 VAL += (*FF)[k+
K]*qq.
data[i+
k] + (*ff)[k+
K]*QQ.
data[i+
k];
919 wno[0].data[
i] +=
val;
920 WNO[0].data[
i] += VAL;
size_t append(const wavearray< DataType_t > &)
virtual void resize(unsigned int)
virtual size_t size() const
void setMatrix(double edge=0., double f=1.)
void _apply_(int n, std::vector< wavearray< double > > &w, std::vector< wavearray< double > > &W)
virtual void rate(double r)
wavearray< double > target
virtual void edge(double s)
TMatrixTSym< double > TMatrixDSym
WDM< double > wdm(nLAYERS, nLAYERS, 6, 10)
std::vector< wavearray< double > > vCROSS
size_t add(WSeries< double > &target, char *name, double fL=0., double fH=0.)
std::vector< double > vectorD
virtual double mean() const
std::vector< wavearray< double > > vEIGEN
std::vector< Wiener > FILTER
std::vector< WSeries< double > > chList
wavearray< double > rank(int nbins=0, double fL=0., double fH=0.)
std::vector< vectorD > filter00
void apply(double threshold=0., char c='a')
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
std::vector< TMatrixDSym > matrix
std::vector< wavearray< int > > chMask
void solve(double th, int nE=0, char c='s')
std::vector< int > channel
wavearray< double > vfreq
std::vector< vectorD > filter90
wavearray< double > getVEIGEN(int n=-1)
std::vector< char * > chName
virtual std::slice getSlice(const double)
std::vector< double > norm
wavearray< double > ts(N)
double fabs(const Complex &x)
wavearray< double > getFILTER(char c='a', int nT=-1, int nW=-1)
void Forward(int n=-1)
param: wavelet - n is number of steps (-1 means full decomposition)
void unmask(int n, double flow=0., double fhigh=0.)
Meyer< double > S(1024, 2)
void mask(int n, double flow=0., double fhigh=0.)
std::vector< wavearray< double > > vrank
WaveDWT< DataType_t > * pWavelet
regression & operator=(const regression &)
wavearray< double > rnoise
void Inverse(int n=-1)
param: n - number of steps (-1 means full reconstruction)
void putLayer(wavearray< DataType_t > &, double n)
param: n - layer number