SHOGUN  3.2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
LogDetEstimator.cpp
Go to the documentation of this file.
1 /*
2  * This program is free software; you can redistribute it and/or modify
3  * it under the terms of the GNU General Public License as published by
4  * the Free Software Foundation; either version 3 of the License, or
5  * (at your option) any later version.
6  *
7  * Written (W) 2013 Soumyajit De
8  */
9 
10 #include <shogun/lib/common.h>
11 #include <shogun/lib/SGVector.h>
12 #include <shogun/lib/SGMatrix.h>
29 
30 namespace shogun
31 {
32 
34  : CSGObject()
35 {
36  init();
37 }
38 
39 #ifdef HAVE_LAPACK
40 #ifdef HAVE_EIGEN3
42  : CSGObject()
43 {
44  init();
45 
46  m_computation_engine=new CSerialComputationEngine();
47  SG_REF(m_computation_engine);
48 
50  new CSparseMatrixOperator<float64_t>(sparse_mat);
51 
52  float64_t accuracy=1E-5;
53 
54  CLanczosEigenSolver* eig_solver=new CLanczosEigenSolver(op);
56 
57  m_operator_log=new CLogRationalApproximationCGM(op,m_computation_engine,
58  eig_solver,linear_solver,accuracy);
59  SG_REF(m_operator_log);
60 
61  #ifdef HAVE_COLPACK
62  m_trace_sampler=new CProbingSampler(op,1,NATURAL,DISTANCE_TWO);
63  #else
64  m_trace_sampler=new CNormalSampler(op->get_dimension());
65  #endif
66 
67  SG_REF(m_trace_sampler);
68 
69  SG_INFO("LogDetEstimator:Using %s, %s with 1E-5 accuracy, %s as default\n",
70  m_computation_engine->get_name(), m_operator_log->get_name(),
71  m_trace_sampler->get_name());
72 }
73 #endif //HAVE_EIGEN3
74 #endif //HAVE_LAPACK
75 
77  COperatorFunction<float64_t>* operator_log,
78  CIndependentComputationEngine* computation_engine)
79  : CSGObject()
80 {
81  init();
82 
83  m_trace_sampler=trace_sampler;
84  SG_REF(m_trace_sampler);
85 
86  m_operator_log=operator_log;
87  SG_REF(m_operator_log);
88 
89  m_computation_engine=computation_engine;
90  SG_REF(m_computation_engine);
91 }
92 
93 void CLogDetEstimator::init()
94 {
95  m_trace_sampler=NULL;
96  m_operator_log=NULL;
97  m_computation_engine=NULL;
98 
99  SG_ADD((CSGObject**)&m_trace_sampler, "trace_sampler",
100  "Trace sampler for the log operator", MS_NOT_AVAILABLE);
101 
102  SG_ADD((CSGObject**)&m_operator_log, "operator_log",
103  "The log operator function", MS_NOT_AVAILABLE);
104 
105  SG_ADD((CSGObject**)&m_computation_engine, "computation_engine",
106  "The computation engine for the jobs", MS_NOT_AVAILABLE);
107 }
108 
110 {
111  SG_UNREF(m_trace_sampler);
112  SG_UNREF(m_operator_log);
113  SG_UNREF(m_computation_engine);
114 }
115 
117 {
118  SG_REF(m_trace_sampler);
119  return m_trace_sampler;
120 }
121 
123 {
124  SG_REF(m_computation_engine);
125  return m_computation_engine;
126 }
127 
129 {
130  SG_REF(m_operator_log);
131  return m_operator_log;
132 }
133 
135 {
136  SG_DEBUG("Entering\n");
137  SG_INFO("Computing %d log-det estimates\n", num_estimates);
138 
139  REQUIRE(m_operator_log, "Operator function is NULL\n");
140  // call the precompute of operator function to compute the prerequisites
141  m_operator_log->precompute();
142 
143  REQUIRE(m_trace_sampler, "Trace sampler is NULL\n");
144  // call the precompute of the sampler
145  m_trace_sampler->precompute();
146 
147  REQUIRE(m_operator_log->get_operator()->get_dimension()\
148  ==m_trace_sampler->get_dimension(),
149  "Mismatch in dimensions of the operator and trace-sampler, %d vs %d!\n",
150  m_operator_log->get_operator()->get_dimension(),
151  m_trace_sampler->get_dimension());
152 
153  // for storing the aggregators that submit_jobs return
154  CDynamicObjectArray* aggregators=new CDynamicObjectArray();
155  index_t num_trace_samples=m_trace_sampler->get_num_samples();
156 
157  for (index_t i=0; i<num_estimates; ++i)
158  {
159  for (index_t j=0; j<num_trace_samples; ++j)
160  {
161  SG_INFO("Computing log-determinant trace sample %d/%d\n", j,
162  num_trace_samples);
163 
164  SG_DEBUG("Creating job for estimate %d, trace sample %d/%d\n", i, j,
165  num_trace_samples);
166  // get the trace sampler vector
167  SGVector<float64_t> s=m_trace_sampler->sample(j);
168  // create jobs with the sample vector and store the aggregator
169  CJobResultAggregator* agg=m_operator_log->submit_jobs(s);
170  aggregators->append_element(agg);
171  SG_UNREF(agg);
172  }
173  }
174 
175  REQUIRE(m_computation_engine, "Computation engine is NULL\n");
176 
177  // wait for all the jobs to be completed
178  SG_INFO("Waiting for jobs to finish\n");
179  m_computation_engine->wait_for_all();
180  SG_INFO("All jobs finished, aggregating results\n");
181 
182  // the samples vector which stores the estimates with averaging
183  SGVector<float64_t> samples(num_estimates);
184  samples.zero();
185 
186  // use the aggregators to find the final result
187  // use the same order as job submission to combine results
188  int32_t num_aggregates=aggregators->get_num_elements();
189  index_t idx_row=0;
190  index_t idx_col=0;
191  for (int32_t i=0; i<num_aggregates; ++i)
192  {
193  // this cast is safe due to above way of building the array
194  CJobResultAggregator* agg=dynamic_cast<CJobResultAggregator*>
195  (aggregators->get_element(i));
196  ASSERT(agg);
197 
198  // call finalize on all the aggregators, cast is safe again
199  agg->finalize();
201  (agg->get_final_result());
202  ASSERT(r);
203 
204  // iterate through indices, group results in the same way as jobs
205  samples[idx_col]+=r->get_result();
206  idx_row++;
207  if (idx_row>=num_trace_samples)
208  {
209  idx_row=0;
210  idx_col++;
211  }
212 
213  SG_UNREF(agg);
214  }
215 
216  // clear all aggregators
217  SG_UNREF(aggregators)
218 
219  SG_INFO("Finished computing %d log-det estimates\n", num_estimates);
220 
221  SG_DEBUG("Leaving\n");
222  return samples;
223 }
224 
226  index_t num_estimates)
227 {
228  SG_DEBUG("Entering...\n")
229 
230  REQUIRE(m_operator_log, "Operator function is NULL\n");
231  // call the precompute of operator function to compute all prerequisites
232  m_operator_log->precompute();
233 
234  REQUIRE(m_trace_sampler, "Trace sampler is NULL\n");
235  // call the precompute of the sampler
236  m_trace_sampler->precompute();
237 
238  // for storing the aggregators that submit_jobs return
239  CDynamicObjectArray aggregators;
240  index_t num_trace_samples=m_trace_sampler->get_num_samples();
241 
242  for (index_t i=0; i<num_estimates; ++i)
243  {
244  for (index_t j=0; j<num_trace_samples; ++j)
245  {
246  // get the trace sampler vector
247  SGVector<float64_t> s=m_trace_sampler->sample(j);
248  // create jobs with the sample vector and store the aggregator
249  CJobResultAggregator* agg=m_operator_log->submit_jobs(s);
250  aggregators.append_element(agg);
251  SG_UNREF(agg);
252  }
253  }
254 
255  REQUIRE(m_computation_engine, "Computation engine is NULL\n");
256  // wait for all the jobs to be completed
257  m_computation_engine->wait_for_all();
258 
259  // the samples matrix which stores the estimates without averaging
260  // dimension: number of trace samples x number of log-det estimates
261  SGMatrix<float64_t> samples(num_trace_samples, num_estimates);
262 
263  // use the aggregators to find the final result
264  int32_t num_aggregates=aggregators.get_num_elements();
265  for (int32_t i=0; i<num_aggregates; ++i)
266  {
267  CJobResultAggregator* agg=dynamic_cast<CJobResultAggregator*>
268  (aggregators.get_element(i));
269  if (!agg)
270  SG_ERROR("Element is not CJobResultAggregator type!\n");
271 
272  // call finalize on all the aggregators
273  agg->finalize();
275  (agg->get_final_result());
276  if (!r)
277  SG_ERROR("Result is not CScalarResult type!\n");
278 
279  // its important that we don't just unref the result here
280  index_t idx_row=i%num_trace_samples;
281  index_t idx_col=i/num_trace_samples;
282  samples(idx_row, idx_col)=r->get_result();
283  SG_UNREF(agg);
284  }
285 
286  // clear all aggregators
287  aggregators.clear_array();
288 
289  SG_DEBUG("Leaving\n")
290  return samples;
291 }
292 
293 }
294 

SHOGUN Machine Learning Toolbox - Documentation