SHOGUN  3.2.1
 全部  命名空间 文件 函数 变量 类型定义 枚举 枚举值 友元 宏定义  
LogDetEstimator.cpp
浏览该文件的文档.
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 
40  : CSGObject()
41 {
42  init();
43 
44  m_computation_engine=new CSerialComputationEngine();
45  SG_REF(m_computation_engine);
46 
48  new CSparseMatrixOperator<float64_t>(sparse_mat);
49 
50  float64_t accuracy=1E-5;
51 
52  CLanczosEigenSolver* eig_solver=new CLanczosEigenSolver(op);
54 
55  m_operator_log=new CLogRationalApproximationCGM(op,m_computation_engine,
56  eig_solver,linear_solver,accuracy);
57  SG_REF(m_operator_log);
58 
59  #ifdef HAVE_COLPACK
60  m_trace_sampler=new CProbingSampler(op,1,NATURAL,DISTANCE_TWO);
61  #else
62  m_trace_sampler=new CNormalSampler(op->get_dimension());
63  #endif
64 
65  SG_REF(m_trace_sampler);
66 
67  SG_INFO("LogDetEstimator:Using %s, %s with 1E-5 accuracy, %s as default\n",
68  m_computation_engine->get_name(), m_operator_log->get_name(),
69  m_trace_sampler->get_name());
70 }
71 
73  COperatorFunction<float64_t>* operator_log,
74  CIndependentComputationEngine* computation_engine)
75  : CSGObject()
76 {
77  init();
78 
79  m_trace_sampler=trace_sampler;
80  SG_REF(m_trace_sampler);
81 
82  m_operator_log=operator_log;
83  SG_REF(m_operator_log);
84 
85  m_computation_engine=computation_engine;
86  SG_REF(m_computation_engine);
87 }
88 
89 void CLogDetEstimator::init()
90 {
91  m_trace_sampler=NULL;
92  m_operator_log=NULL;
93  m_computation_engine=NULL;
94 
95  SG_ADD((CSGObject**)&m_trace_sampler, "trace_sampler",
96  "Trace sampler for the log operator", MS_NOT_AVAILABLE);
97 
98  SG_ADD((CSGObject**)&m_operator_log, "operator_log",
99  "The log operator function", MS_NOT_AVAILABLE);
100 
101  SG_ADD((CSGObject**)&m_computation_engine, "computation_engine",
102  "The computation engine for the jobs", MS_NOT_AVAILABLE);
103 }
104 
106 {
107  SG_UNREF(m_trace_sampler);
108  SG_UNREF(m_operator_log);
109  SG_UNREF(m_computation_engine);
110 }
111 
113 {
114  SG_REF(m_trace_sampler);
115  return m_trace_sampler;
116 }
117 
119 {
120  SG_REF(m_computation_engine);
121  return m_computation_engine;
122 }
123 
125 {
126  SG_REF(m_operator_log);
127  return m_operator_log;
128 }
129 
131 {
132  SG_DEBUG("Entering\n");
133  SG_INFO("Computing %d log-det estimates\n", num_estimates);
134 
135  REQUIRE(m_operator_log, "Operator function is NULL\n");
136  // call the precompute of operator function to compute the prerequisites
137  m_operator_log->precompute();
138 
139  REQUIRE(m_trace_sampler, "Trace sampler is NULL\n");
140  // call the precompute of the sampler
141  m_trace_sampler->precompute();
142 
143  REQUIRE(m_operator_log->get_operator()->get_dimension()\
144  ==m_trace_sampler->get_dimension(),
145  "Mismatch in dimensions of the operator and trace-sampler, %d vs %d!\n",
146  m_operator_log->get_operator()->get_dimension(),
147  m_trace_sampler->get_dimension());
148 
149  // for storing the aggregators that submit_jobs return
150  CDynamicObjectArray* aggregators=new CDynamicObjectArray();
151  index_t num_trace_samples=m_trace_sampler->get_num_samples();
152 
153  for (index_t i=0; i<num_estimates; ++i)
154  {
155  for (index_t j=0; j<num_trace_samples; ++j)
156  {
157  SG_INFO("Computing log-determinant trace sample %d/%d\n", j,
158  num_trace_samples);
159 
160  SG_DEBUG("Creating job for estimate %d, trace sample %d/%d\n", i, j,
161  num_trace_samples);
162  // get the trace sampler vector
163  SGVector<float64_t> s=m_trace_sampler->sample(j);
164  // create jobs with the sample vector and store the aggregator
165  CJobResultAggregator* agg=m_operator_log->submit_jobs(s);
166  aggregators->append_element(agg);
167  SG_UNREF(agg);
168  }
169  }
170 
171  REQUIRE(m_computation_engine, "Computation engine is NULL\n");
172 
173  // wait for all the jobs to be completed
174  SG_INFO("Waiting for jobs to finish\n");
175  m_computation_engine->wait_for_all();
176  SG_INFO("All jobs finished, aggregating results\n");
177 
178  // the samples vector which stores the estimates with averaging
179  SGVector<float64_t> samples(num_estimates);
180  samples.zero();
181 
182  // use the aggregators to find the final result
183  // use the same order as job submission to combine results
184  int32_t num_aggregates=aggregators->get_num_elements();
185  index_t idx_row=0;
186  index_t idx_col=0;
187  for (int32_t i=0; i<num_aggregates; ++i)
188  {
189  // this cast is safe due to above way of building the array
190  CJobResultAggregator* agg=dynamic_cast<CJobResultAggregator*>
191  (aggregators->get_element(i));
192  ASSERT(agg);
193 
194  // call finalize on all the aggregators, cast is safe again
195  agg->finalize();
197  (agg->get_final_result());
198  ASSERT(r);
199 
200  // iterate through indices, group results in the same way as jobs
201  samples[idx_col]+=r->get_result();
202  idx_row++;
203  if (idx_row>=num_trace_samples)
204  {
205  idx_row=0;
206  idx_col++;
207  }
208 
209  SG_UNREF(agg);
210  }
211 
212  // clear all aggregators
213  SG_UNREF(aggregators)
214 
215  SG_INFO("Finished computing %d log-det estimates\n", num_estimates);
216 
217  SG_DEBUG("Leaving\n");
218  return samples;
219 }
220 
222  index_t num_estimates)
223 {
224  SG_DEBUG("Entering...\n")
225 
226  REQUIRE(m_operator_log, "Operator function is NULL\n");
227  // call the precompute of operator function to compute all prerequisites
228  m_operator_log->precompute();
229 
230  REQUIRE(m_trace_sampler, "Trace sampler is NULL\n");
231  // call the precompute of the sampler
232  m_trace_sampler->precompute();
233 
234  // for storing the aggregators that submit_jobs return
235  CDynamicObjectArray aggregators;
236  index_t num_trace_samples=m_trace_sampler->get_num_samples();
237 
238  for (index_t i=0; i<num_estimates; ++i)
239  {
240  for (index_t j=0; j<num_trace_samples; ++j)
241  {
242  // get the trace sampler vector
243  SGVector<float64_t> s=m_trace_sampler->sample(j);
244  // create jobs with the sample vector and store the aggregator
245  CJobResultAggregator* agg=m_operator_log->submit_jobs(s);
246  aggregators.append_element(agg);
247  SG_UNREF(agg);
248  }
249  }
250 
251  REQUIRE(m_computation_engine, "Computation engine is NULL\n");
252  // wait for all the jobs to be completed
253  m_computation_engine->wait_for_all();
254 
255  // the samples matrix which stores the estimates without averaging
256  // dimension: number of trace samples x number of log-det estimates
257  SGMatrix<float64_t> samples(num_trace_samples, num_estimates);
258 
259  // use the aggregators to find the final result
260  int32_t num_aggregates=aggregators.get_num_elements();
261  for (int32_t i=0; i<num_aggregates; ++i)
262  {
263  CJobResultAggregator* agg=dynamic_cast<CJobResultAggregator*>
264  (aggregators.get_element(i));
265  if (!agg)
266  SG_ERROR("Element is not CJobResultAggregator type!\n");
267 
268  // call finalize on all the aggregators
269  agg->finalize();
271  (agg->get_final_result());
272  if (!r)
273  SG_ERROR("Result is not CScalarResult type!\n");
274 
275  // its important that we don't just unref the result here
276  index_t idx_row=i%num_trace_samples;
277  index_t idx_col=i/num_trace_samples;
278  samples(idx_row, idx_col)=r->get_result();
279  SG_UNREF(agg);
280  }
281 
282  // clear all aggregators
283  aggregators.clear_array();
284 
285  SG_DEBUG("Leaving\n")
286  return samples;
287 }
288 
289 }
290 
Base class that stores the result of an independent job when the result is a scalar.
template class SGSparseMatrix
#define SG_INFO(...)
Definition: SGIO.h:120
virtual const index_t get_dimension() const
Definition: TraceSampler.h:77
CTraceSampler * get_trace_sampler(void) const
CLinearOperator< T > * get_operator() const
SGVector< float64_t > sample(index_t num_estimates)
int32_t index_t
Definition: common.h:60
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...
#define SG_UNREF(x)
Definition: SGRefObject.h:35
virtual SGVector< float64_t > sample(index_t idx) const =0
#define SG_ERROR(...)
Definition: SGIO.h:131
#define REQUIRE(x,...)
Definition: SGIO.h:208
virtual const index_t get_num_samples() const
Definition: TraceSampler.h:71
SGMatrix< float64_t > sample_without_averaging(index_t num_estimates)
virtual const char * get_name() const
CJobResult * get_final_result() const
virtual void precompute()=0
#define ASSERT(x)
Definition: SGIO.h:203
Class SGObject is the base class of all shogun objects.
Definition: SGObject.h:102
Class that computes eigenvalues of a real valued, self-adjoint linear operator using Lanczos algorith...
double float64_t
Definition: common.h:48
virtual const char * get_name() const
Definition: TraceSampler.h:83
#define SG_REF(x)
Definition: SGRefObject.h:34
CIndependentComputationEngine * get_computation_engine(void) const
virtual CJobResultAggregator * submit_jobs(SGVector< T > sample)=0
Dynamic array class for CSGObject pointers that creates an array that can be used like a list or an a...
virtual void precompute()=0
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
const index_t get_dimension() const
#define SG_DEBUG(...)
Definition: SGIO.h:109
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:66
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:71
bool append_element(CSGObject *e)

SHOGUN 机器学习工具包 - 项目文档