3 #pragma GCC system_header
10 #include "TObjArray.h"
11 #include "TObjString.h"
35 cout <<
"-----> CWB_Plugin_Linear_Bilinear_Regression_Rank.C : " << ifo.Data() << endl;
42 bool SAVE_TREE =
false;
100 IMPORT(
int,RESAMPLING_INDEX)
109 if(FRLIST_WITNESS==
"") {cout <<
"Error : witness frames list not defined" << endl;gSystem->Exit(1);}
110 if(CHLIST_LINEAR==
"") {cout <<
"Error : linear witness channels list not defined" << endl;gSystem->Exit(1);}
111 if(CHLIST_WITNESS==
"") {cout <<
"Error : witness channels list not defined" << endl;gSystem->Exit(1);}
112 if(CHNAME_LINEAR==
"") {cout <<
"Error : main linear witness channel list not defined" << endl;gSystem->Exit(1);}
114 IMPORT(
double,LAYER_WIDTH)
121 IMPORT(
double,L_APPLY_THRESHOLD)
122 IMPORT(
double,L_SOLVE_THRESHOLD)
123 IMPORT(
double,L_SOLVE_NEIGEN_PER_LAYER)
124 IMPORT(
char,L_SOLVE_REGULATOR)
127 IMPORT(
double,B_APPLY_THRESHOLD)
128 IMPORT(
double,B_SOLVE_THRESHOLD)
129 IMPORT(
double,B_SOLVE_NEIGEN_PER_LAYER)
130 IMPORT(
char,B_SOLVE_REGULATOR)
132 REGRESSION_RANK.ToUpper();
133 if(REGRESSION_RANK!=
"LINEAR" && REGRESSION_RANK!=
"BILINEAR") {
134 cout <<
"Error : REGRESSION_RANK not defined [LINEAR/BILINEAR]" << endl;
140 int rank[4] = {50,30,10,1};
143 std::vector<std::string> CHN;
144 std::vector<std::string>* pCHN=&CHN;
149 for(
int n=0;
n<cfg->
nIFO;
n++)
if(ifo==cfg->
ifo[
n]) {xIFO=
n;
break;}
151 int runID = net->
nRun;
175 while (((SR % 2) == 0) && SR > 2*
CUT_LOW_FREQ) {SR /= 2;mm++;}
180 cout <<
"Set Frequency Range [0:" << SR/2 <<
"] = 0" << endl;
193 if(scratch>cfg->
segEdge+0.001) {
195 cout <<
"Regression Plugin : Error - filter scratch must be <= cwb scratch!!!" << endl;
196 cout <<
"filter scratch : " << scratch <<
" sec" << endl;
197 cout <<
"cwb scratch : " << cfg->
segEdge <<
" sec" << endl;
207 cout <<
"CWB_PLUGIN_DATA_MDC : Apply Linear Regression ..." << endl;
209 frame frl(FRLIST_WITNESS,
"",
"README",
true);
214 if(REGRESSION_RANK==
"LINEAR") {
215 cout <<
"CWB_PLUGIN_DATA_MDC : Apply Linear Regression Rank ..." << endl;
219 ifl.open(CHLIST_WITNESS.Data(),
ios::in);
220 if (!ifl.good()) {cout <<
"Error Opening File : " << CHLIST_WITNESS.Data() << endl;gSystem->Exit(1);}
226 sprintf(ofname,
"%s/linear_regression_rank_%s_%s.txt",cfg->
output_dir,ifo.Data(),flabel);
227 cout <<
"write results to " << ofname << endl;
231 if (!out.good()) {cout <<
"Error Opening File : " << ofname << endl;gSystem->Exit(1);}
237 if (!ifl.good())
break;
238 if(witness[0]==
'#')
continue;
243 size = rr.
add(wT,
"target");
244 if(size==0) {cout <<
"Regression Plugin - empty target channel" << endl;gSystem->Exit(1);}
247 size = rr.
add(wl,witness);
251 rr.
solve(L_SOLVE_THRESHOLD,L_SOLVE_NEIGEN_PER_LAYER,L_SOLVE_REGULATOR);
252 rr.
apply(L_APPLY_THRESHOLD);
254 for(
int n=0;
n<4;
n++) {
255 e[
n] = rr.
rank(rank[
n],0);
259 for(
int n=0;
n<4;
n++) {
263 CHN.push_back(witness);
264 sprintf(fstring,
"%3d %35s %10.5f %9.4f %8.3f %7.2f",
265 nCH,witness,CHC[0][nCH],CHC[1][nCH],CHC[2][nCH],CHC[3][nCH]);
266 cout << fstring << endl;
267 if(SAVE_ASCII) out << fstring << endl;out.flush();
272 if(SAVE_ASCII) out.close();
277 ifl.open(CHLIST_LINEAR.Data(),
ios::in);
278 if (!ifl.good()) {cout <<
"Error Opening File : " << CHLIST_LINEAR.Data() << endl;gSystem->Exit(1);}
283 size = rr.
add(wT,
"target");
284 if(size==0) {cout <<
"Regression Plugin - empty target channel" << endl;gSystem->Exit(1);}
289 if (!ifl.good())
break;
290 if(linear[0]==
'#')
continue;
291 cout <<
"read linear channel : \t" << linear << endl;
301 rr.
solve(L_SOLVE_THRESHOLD,L_SOLVE_NEIGEN_PER_LAYER,L_SOLVE_REGULATOR);
302 rr.
apply(L_APPLY_THRESHOLD);
310 if(REGRESSION_RANK==
"BILINEAR") {
312 cout <<
"CWB_PLUGIN_DATA_MDC : Apply Bilinear Regression Rank ..." << endl;
315 sprintf(ofname,
"%s/bilinear_regression_rank_%s_%s.txt",cfg->
output_dir,ifo.Data(),flabel);
316 cout <<
"write results to " << ofname << endl;
320 if (!out.good()) {cout <<
"Error Opening File : " << ofname << endl;gSystem->Exit(1);}
328 frame frb(FRLIST_WITNESS,
"",
"README",
true);
332 iwb.open(CHLIST_WITNESS.Data(),
ios::in);
333 if (!iwb.good()) {cout <<
"Error Opening File : " << CHLIST_WITNESS.Data() << endl;gSystem->Exit(2);}
343 if (!iwb.good())
break;
345 if(witness[0]==
'#')
continue;
351 size=rr.
add(wT,
"target");
352 if(size==0) {cout <<
"Regression Plugin - empty target channel" << endl;gSystem->Exit(1);}
363 if(size==0) {cout <<
"Regression Plugin - empty main linear channel" << endl;gSystem->Exit(1);}
365 if(rr.
add(wb,witness,WFLOW,WFHIGH)) {
366 rr.
add(1,2,
"biwitness");
371 rr.
solve(B_SOLVE_THRESHOLD,B_SOLVE_NEIGEN_PER_LAYER,B_SOLVE_REGULATOR);
372 rr.
apply(B_APPLY_THRESHOLD);
374 for(
int n=0;
n<4;
n++) {
375 e[
n] = rr.
rank(rank[
n],0);
379 for(
int n=0;
n<4;
n++) {
383 CHN.push_back(witness);
385 sprintf(fstring,
"%3d %35s %10.5f %9.4f %8.3f %7.2f",
386 nCH,witness,CHC[0][nCH],CHC[1][nCH],CHC[2][nCH],CHC[3][nCH]);
387 cout << fstring << endl;
388 if(SAVE_ASCII) out << fstring << endl;out.flush();
391 if(SAVE_ASCII) out.close();
401 if(REGRESSION_RANK==
"LINEAR") {
402 sprintf(outFile,
"%s/linear_regression_rank_%s_%s.root",cfg->
output_dir,ifo.Data(),flabel);
404 sprintf(outFile,
"%s/bilinear_regression_rank_%s_%s.root",cfg->
output_dir,ifo.Data(),flabel);
406 cout <<
"write results to " << outFile << endl;
408 TFile
ofile = TFile(outFile,
"RECREATE");
409 cfg->Write(
"config");
412 TTree rrtree(
"regression",
"regression");
413 double igps=x->
start();
416 char sifo[256];
sprintf(sifo,
"%s",ifo.Data());
417 rrtree.Branch(
"run",&runID,
"run/I");
418 rrtree.Branch(
"gps",&igps,
"gps/D");
419 rrtree.Branch(
"start",&start,
"start/D");
420 rrtree.Branch(
"stop",&stop,
"stop/D");
421 rrtree.Branch(
"sifo",sifo,
"sifo/C");
422 rrtree.Branch(
"ifo",&xIFO,
"ifo/I");
423 rrtree.Branch(
"nch",&nCH,
"nch/I");
424 rrtree.Branch(
"chc50",CHC[0],
"chc[nch]/D");
425 rrtree.Branch(
"chc30",CHC[1],
"chc[nch]/D");
426 rrtree.Branch(
"chc10",CHC[2],
"chc[nch]/D");
427 rrtree.Branch(
"chc01",CHC[3],
"chc[nch]/D");
428 rrtree.Branch(
"chn",&pCHN);
436 if(ifo==cfg->
ifo[xIFO]) {
437 cout <<
"remove temporary file ..." << endl;
439 jname.ReplaceAll(
":/",
"");
440 cout << jname.Data() << endl;
441 gSystem->Exec(
TString(
"rm "+jname).Data());
442 cout <<
"end job" << endl;
double L_SOLVE_NEIGEN_PER_LAYER
virtual size_t size() const
void setMatrix(double edge=0., double f=1.)
virtual void rate(double r)
cout<< "skymap size : "<< L<< endl;for(int l=0;l< L;l++) sm.set(l, l);sm > const_cast< char * >("skymap.dat")
double B_SOLVE_NEIGEN_PER_LAYER
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\t layers : "<< nLAYERS<< "\t dF(hz) : "<< dF<< "\t dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1)*itime+ifreq;double time=itime *dT;double freq=(ifreq >0)?ifreq *dF:dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
size_t add(WSeries< double > &target, char *name, double fL=0., double fH=0.)
virtual void start(double s)
network ** net
NOISE_MDC_SIMULATION.
void setChName(TString chName)
wavearray< double > rank(int nbins=0, double fL=0., double fH=0.)
#define IMPORT(TYPE, VAR)
void apply(double threshold=0., char c='a')
int getLayer(wavearray< DataType_t > &w, double n)
param: n - layer number
void CWB_Plugin(TFile *jfile, CWB::config *cfg, network *net, WSeries< double > *x, TString ifo, int type)
COHERENCE.
void solve(double th, int nE=0, char c='s')
void setSRIndex(int srIndex)
wavearray< double > getClean()
virtual void stop(double s)
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.)
void mask(int n, double flow=0., double fhigh=0.)
sprintf(tfres,"(1/%g)x(%g) (sec)x(Hz)", 2 *df, df)
virtual void resize(unsigned int)
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