test_ress.cxx
Go to the documentation of this file.
3 
4 #include <Eigen/Dense>
5 using namespace Eigen;
6 
7 #include <iostream>
8 using namespace std;
9 
10 void test_model(WireCell::LinearModel& m, MatrixXd& G, VectorXd& W);
11 void test_lasso(WireCell::LassoModel& m, MatrixXd& G, VectorXd& W);
12 void print_results(WireCell::LinearModel& m, VectorXd& C);
13 
14 int main(int argc, char* argv[])
15 {
16  // std::srand((unsigned int) time(0));
17 
18  const int N_CELL = 40;
19  const int N_ZERO = 30;
20  const int N_WIRE = int(N_CELL * 0.7);
21 
22  // initialize C vector: NCELL cells with NZERO zeros. (true charge in each cell)
23  VectorXd C = VectorXd::Random(N_CELL)*50 + VectorXd::Constant(N_CELL, 150);
24  VectorXd r = N_CELL / 2 * (VectorXd::Random(N_ZERO)+VectorXd::Constant(N_ZERO, 1));
25  for (int i=0; i<N_ZERO; i++) {
26  C( int(r(i)) ) = 0;
27  }
28 
29  // initialize G matrix: N_WIRE rows and N_CELL columns. (geometry matrix)
30  MatrixXd G = MatrixXd::Zero(N_WIRE, N_CELL);
31  for (int i=0; i<N_WIRE; i++) {
32  VectorXd t = VectorXd::Random(N_CELL);
33  for (int j=0; j<N_CELL; j++) {
34  G(i, j) = int(t(j)+1);
35  }
36  }
37 
38  // W vector is the measured charge on wires.
39  VectorXd W = G * C;
40 
41  cout << "geometry matrix:" << endl;
42  cout << G << endl << endl;
43 
44  // cout << "measured charge on each wire:" << endl;
45  // cout << W.transpose() << endl << endl;
46 
47  cout << "true charge of each cell:" << endl;
48  cout << C.transpose() << endl << endl;
49  double lambda = 0.01;
50 
51  // WireCell::ElasticNetModel m(lambda, 0.95, 100000, 1e-4);
52  // test_model(m, G, W);
53  // print_results(m, C);
54 
55  WireCell::LassoModel m2(lambda, 100000, 1e-3);
56  test_lasso(m2, G, W);
57  print_results(m2, C);
58 
59  return 0;
60 }
61 
62 void test_model(WireCell::LinearModel& m, MatrixXd& G, VectorXd& W)
63 {
64  m.SetData(G, W);
65  m.Fit();
66 
67 }
68 
69 void test_lasso(WireCell::LassoModel& m, MatrixXd& G, VectorXd& W)
70 {
71  m.SetData(G, W);
72 
73  // one can set the weight of each cell's regularization.
74  // m.SetLambdaWeight(0, 10.);
75  // m.SetLambdaWeight(2, 10.);
76  // m.SetLambdaWeight(6, 10.);
77  m.Fit();
78 
79  cout << "chi2: " << m.name << ": base="
80  << m.chi2_base() << ", l1=" << m.chi2_l1()<< endl << endl;
81 }
82 
84 {
85  VectorXd beta = m.Getbeta();
86 
87  cout << "fitted charge of each cell: " << m.name << endl;
88  cout << beta.transpose() << endl << endl;
89 
90  // cout << "predicted charge on each wire: Lasso" << endl;
91  // cout << m.Predict().transpose() << endl << endl;
92 
93  cout << "average residual charge difference per wire: " << m.name << ": "
94  << m.MeanResidual() << endl;
95 
96 
97  int nbeta = beta.size();
98 
99  int n_zero_true = 0;
100  int n_zero_beta = 0;
101  int n_zero_correct = 0;
102 
103  for (int i=0; i<nbeta; i++) {
104  if (fabs(C(i))<0.1) n_zero_true++;
105  if (fabs(beta(i))<5) n_zero_beta++;
106  if (fabs(C(i))<0.1 && (fabs(C(i) - beta(i)) < 10)) n_zero_correct++;
107  }
108  cout << "true zeros: " << n_zero_true << endl;
109  cout << "fitted zeros: " << n_zero_beta << endl;
110  cout << "correct fitted zeros: " << n_zero_correct << endl;
111 
112  cout << endl;
113 }
static const double m
Definition: Units.h:79
STL namespace.
void print_results(WireCell::LinearModel &m, VectorXd &C)
Definition: test_ress.cxx:83
const double e
Eigen::VectorXd & Getbeta()
Definition: LinearModel.h:16
void test_model(WireCell::LinearModel &m, MatrixXd &G, VectorXd &W)
Definition: test_ress.cxx:62
virtual void SetData(Eigen::MatrixXd X, Eigen::VectorXd y)
Definition: LinearModel.h:18
int main(int argc, char *argv[])
Definition: test_ress.cxx:14
static const double m2
Definition: Units.h:80
void test_lasso(WireCell::LassoModel &m, MatrixXd &G, VectorXd &W)
Definition: test_ress.cxx:69
virtual void Fit()
Definition: LinearModel.h:23
QTextStream & endl(QTextStream &s)