SHOGUN  5.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
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
41  : CSGObject()
42 {
43  init();
44 
45  m_computation_engine=new CSerialComputationEngine();
46  SG_REF(m_computation_engine);
47 
49  new CSparseMatrixOperator<float64_t>(sparse_mat);
50 
51  float64_t accuracy=1E-5;
52 
53  CLanczosEigenSolver* eig_solver=new CLanczosEigenSolver(op);
55 
56  m_operator_log=new CLogRationalApproximationCGM(op,m_computation_engine,
57  eig_solver,linear_solver,accuracy);
58  SG_REF(m_operator_log);
59 
60  #ifdef HAVE_COLPACK
61  m_trace_sampler=new CProbingSampler(op,1,NATURAL,DISTANCE_TWO);
62  #else
63  m_trace_sampler=new CNormalSampler(op->get_dimension());
64  #endif
65 
66  SG_REF(m_trace_sampler);
67 
68  SG_INFO("LogDetEstimator:Using %s, %s with 1E-5 accuracy, %s as default\n",
69  m_computation_engine->get_name(), m_operator_log->get_name(),
70  m_trace_sampler->get_name());
71 }
72 #endif //HAVE_LAPACK
73 
75  COperatorFunction<float64_t>* operator_log,
76  CIndependentComputationEngine* computation_engine)
77  : CSGObject()
78 {
79  init();
80 
81  m_trace_sampler=trace_sampler;
82  SG_REF(m_trace_sampler);
83 
84  m_operator_log=operator_log;
85  SG_REF(m_operator_log);
86 
87  m_computation_engine=computation_engine;
88  SG_REF(m_computation_engine);
89 }
90 
91 void CLogDetEstimator::init()
92 {
93  m_trace_sampler=NULL;
94  m_operator_log=NULL;
95  m_computation_engine=NULL;
96 
97  SG_ADD((CSGObject**)&m_trace_sampler, "trace_sampler",
98  "Trace sampler for the log operator", MS_NOT_AVAILABLE);
99 
100  SG_ADD((CSGObject**)&m_operator_log, "operator_log",
101  "The log operator function", MS_NOT_AVAILABLE);
102 
103  SG_ADD((CSGObject**)&m_computation_engine, "computation_engine",
104  "The computation engine for the jobs", MS_NOT_AVAILABLE);
105 }
106 
108 {
109  SG_UNREF(m_trace_sampler);
110  SG_UNREF(m_operator_log);
111  SG_UNREF(m_computation_engine);
112 }
113 
115 {
116  SG_REF(m_trace_sampler);
117  return m_trace_sampler;
118 }
119 
121 {
122  SG_REF(m_computation_engine);
123  return m_computation_engine;
124 }
125 
127 {
128  SG_REF(m_operator_log);
129  return m_operator_log;
130 }
131 
133 {
134  SG_DEBUG("Entering\n");
135  SG_INFO("Computing %d log-det estimates\n", num_estimates);
136 
137  REQUIRE(m_operator_log, "Operator function is NULL\n");
138  // call the precompute of operator function to compute the prerequisites
139  m_operator_log->precompute();
140 
141  REQUIRE(m_trace_sampler, "Trace sampler is NULL\n");
142  // call the precompute of the sampler
143  m_trace_sampler->precompute();
144 
145  REQUIRE(m_operator_log->get_operator()->get_dimension()\
146  ==m_trace_sampler->get_dimension(),
147  "Mismatch in dimensions of the operator and trace-sampler, %d vs %d!\n",
148  m_operator_log->get_operator()->get_dimension(),
149  m_trace_sampler->get_dimension());
150 
151  // for storing the aggregators that submit_jobs return
152  CDynamicObjectArray* aggregators=new CDynamicObjectArray();
153  index_t num_trace_samples=m_trace_sampler->get_num_samples();
154 
155  for (index_t i=0; i<num_estimates; ++i)
156  {
157  for (index_t j=0; j<num_trace_samples; ++j)
158  {
159  SG_INFO("Computing log-determinant trace sample %d/%d\n", j,
160  num_trace_samples);
161 
162  SG_DEBUG("Creating job for estimate %d, trace sample %d/%d\n", i, j,
163  num_trace_samples);
164  // get the trace sampler vector
165  SGVector<float64_t> s=m_trace_sampler->sample(j);
166  // create jobs with the sample vector and store the aggregator
167  CJobResultAggregator* agg=m_operator_log->submit_jobs(s);
168  aggregators->append_element(agg);
169  SG_UNREF(agg);
170  }
171  }
172 
173  REQUIRE(m_computation_engine, "Computation engine is NULL\n");
174 
175  // wait for all the jobs to be completed
176  SG_INFO("Waiting for jobs to finish\n");
177  m_computation_engine->wait_for_all();
178  SG_INFO("All jobs finished, aggregating results\n");
179 
180  // the samples vector which stores the estimates with averaging
181  SGVector<float64_t> samples(num_estimates);
182  samples.zero();
183 
184  // use the aggregators to find the final result
185  // use the same order as job submission to combine results
186  int32_t num_aggregates=aggregators->get_num_elements();
187  index_t idx_row=0;
188  index_t idx_col=0;
189  for (int32_t i=0; i<num_aggregates; ++i)
190  {
191  // this cast is safe due to above way of building the array
192  CJobResultAggregator* agg=dynamic_cast<CJobResultAggregator*>
193  (aggregators->get_element(i));
194  ASSERT(agg);
195 
196  // call finalize on all the aggregators, cast is safe again
197  agg->finalize();
199  (agg->get_final_result());
200  ASSERT(r);
201 
202  // iterate through indices, group results in the same way as jobs
203  samples[idx_col]+=r->get_result();
204  idx_row++;
205  if (idx_row>=num_trace_samples)
206  {
207  idx_row=0;
208  idx_col++;
209  }
210 
211  SG_UNREF(agg);
212  }
213 
214  // clear all aggregators
215  SG_UNREF(aggregators)
216 
217  SG_INFO("Finished computing %d log-det estimates\n", num_estimates);
218 
219  SG_DEBUG("Leaving\n");
220  return samples;
221 }
222 
224  index_t num_estimates)
225 {
226  SG_DEBUG("Entering...\n")
227 
228  REQUIRE(m_operator_log, "Operator function is NULL\n");
229  // call the precompute of operator function to compute all prerequisites
230  m_operator_log->precompute();
231 
232  REQUIRE(m_trace_sampler, "Trace sampler is NULL\n");
233  // call the precompute of the sampler
234  m_trace_sampler->precompute();
235 
236  // for storing the aggregators that submit_jobs return
237  CDynamicObjectArray aggregators;
238  index_t num_trace_samples=m_trace_sampler->get_num_samples();
239 
240  for (index_t i=0; i<num_estimates; ++i)
241  {
242  for (index_t j=0; j<num_trace_samples; ++j)
243  {
244  // get the trace sampler vector
245  SGVector<float64_t> s=m_trace_sampler->sample(j);
246  // create jobs with the sample vector and store the aggregator
247  CJobResultAggregator* agg=m_operator_log->submit_jobs(s);
248  aggregators.append_element(agg);
249  SG_UNREF(agg);
250  }
251  }
252 
253  REQUIRE(m_computation_engine, "Computation engine is NULL\n");
254  // wait for all the jobs to be completed
255  m_computation_engine->wait_for_all();
256 
257  // the samples matrix which stores the estimates without averaging
258  // dimension: number of trace samples x number of log-det estimates
259  SGMatrix<float64_t> samples(num_trace_samples, num_estimates);
260 
261  // use the aggregators to find the final result
262  int32_t num_aggregates=aggregators.get_num_elements();
263  for (int32_t i=0; i<num_aggregates; ++i)
264  {
265  CJobResultAggregator* agg=dynamic_cast<CJobResultAggregator*>
266  (aggregators.get_element(i));
267  if (!agg)
268  SG_ERROR("Element is not CJobResultAggregator type!\n");
269 
270  // call finalize on all the aggregators
271  agg->finalize();
273  (agg->get_final_result());
274  if (!r)
275  SG_ERROR("Result is not CScalarResult type!\n");
276 
277  // its important that we don't just unref the result here
278  index_t idx_row=i%num_trace_samples;
279  index_t idx_col=i/num_trace_samples;
280  samples(idx_row, idx_col)=r->get_result();
281  SG_UNREF(agg);
282  }
283 
284  // clear all aggregators
285  aggregators.clear_array();
286 
287  SG_DEBUG("Leaving\n")
288  return samples;
289 }
290 
291 }
292 
Base class that stores the result of an independent job when the result is a scalar.
Definition: ScalarResult.h:24
#define SG_INFO(...)
Definition: SGIO.h:118
virtual const index_t get_dimension() const
Definition: TraceSampler.h:77
CTraceSampler * get_trace_sampler(void) const
const index_t get_dimension() const
SGVector< float64_t > sample(index_t num_estimates)
int32_t index_t
Definition: common.h:62
Class that computes multiple independent instances of computation jobs sequentially.
Implementaion of rational approximation of a operator-function times vector where the operator functi...
class that uses conjugate gradient method for solving a shifted linear system family where the linear...
virtual SGVector< float64_t > sample(index_t idx) const =0
virtual void precompute()=0
#define SG_ERROR(...)
Definition: SGIO.h:129
#define REQUIRE(x,...)
Definition: SGIO.h:206
virtual const index_t get_num_samples() const
Definition: TraceSampler.h:71
SGMatrix< float64_t > sample_without_averaging(index_t num_estimates)
#define SG_REF(x)
Definition: SGObject.h:54
CJobResult * get_final_result() const
virtual void precompute()=0
#define ASSERT(x)
Definition: SGIO.h:201
Class SGObject is the base class of all shogun objects.
Definition: SGObject.h:115
Class that computes eigenvalues of a real valued, self-adjoint linear operator using Lanczos algorith...
double float64_t
Definition: common.h:50
virtual const char * get_name() const
Definition: TraceSampler.h:83
CLinearOperator< T > * get_operator() const
CIndependentComputationEngine * get_computation_engine(void) const
virtual const char * get_name() const
Dynamic array class for CSGObject pointers that creates an array that can be used like a list or an a...
Abstract base class that provides an interface for computing an aggeregation of the job results of in...
COperatorFunction< float64_t > * get_operator_function(void) const
#define SG_UNREF(x)
Definition: SGObject.h:55
#define SG_DEBUG(...)
Definition: SGIO.h:107
all of classes and functions are contained in the shogun namespace
Definition: class_list.h:18
Class that represents a sparse-matrix linear operator. It computes matrix-vector product in its appl...
Class that provides a sample method for Gaussian samples.
Definition: NormalSampler.h:22
Abstract base class for solving multiple independent instances of CIndependentJob. It has one method, submit_job, which may add the job to an internal queue and might block if there is yet not space in the queue. After jobs are submitted, it might not yet be ready. wait_for_all waits until all jobs are completed, which *must be* called to guarantee that all jobs are finished.
const T get_result() const
Definition: ScalarResult.h:67
CSGObject * get_element(int32_t index) const
Abstract template base class that provides an interface for sampling the trace of a linear operator u...
Definition: TraceSampler.h:23
#define SG_ADD(...)
Definition: SGObject.h:84
virtual CJobResultAggregator * submit_jobs(SGVector< T > sample)=0
bool append_element(CSGObject *e)

SHOGUN Machine Learning Toolbox - Documentation