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 

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