SHOGUN  4.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules 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 
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:51
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:112
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:52
#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:81
virtual CJobResultAggregator * submit_jobs(SGVector< T > sample)=0
bool append_element(CSGObject *e)

SHOGUN Machine Learning Toolbox - Documentation