SHOGUN  v3.2.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
LogDetEstimator.cpp
Go to the documentation of this file.
00001 /*
00002  * This program is free software; you can redistribute it and/or modify
00003  * it under the terms of the GNU General Public License as published by
00004  * the Free Software Foundation; either version 3 of the License, or
00005  * (at your option) any later version.
00006  *
00007  * Written (W) 2013 Soumyajit De
00008  */
00009 
00010 #include <shogun/lib/common.h>
00011 #include <shogun/lib/SGVector.h>
00012 #include <shogun/lib/SGMatrix.h>
00013 #include <shogun/lib/DynamicObjectArray.h>
00014 #include <shogun/lib/computation/engine/IndependentComputationEngine.h>
00015 #include <shogun/lib/computation/jobresult/ScalarResult.h>
00016 #include <shogun/lib/computation/aggregator/JobResultAggregator.h>
00017 #include <shogun/mathematics/linalg/ratapprox/tracesampler/TraceSampler.h>
00018 #include <shogun/mathematics/linalg/ratapprox/opfunc/OperatorFunction.h>
00019 #include <shogun/mathematics/linalg/ratapprox/logdet/LogDetEstimator.h>
00020 
00021 namespace shogun
00022 {
00023 
00024 CLogDetEstimator::CLogDetEstimator()
00025     : CSGObject()
00026 {
00027     init();
00028 }
00029 
00030 CLogDetEstimator::CLogDetEstimator(CTraceSampler* trace_sampler,
00031     COperatorFunction<float64_t>* operator_log,
00032     CIndependentComputationEngine* computation_engine)
00033     : CSGObject()
00034 {
00035     init();
00036 
00037     m_trace_sampler=trace_sampler;
00038     SG_REF(m_trace_sampler);
00039 
00040     m_operator_log=operator_log;
00041     SG_REF(m_operator_log);
00042 
00043     m_computation_engine=computation_engine;
00044     SG_REF(m_computation_engine);
00045 }
00046 
00047 void CLogDetEstimator::init()
00048 {
00049     m_trace_sampler=NULL;
00050     m_operator_log=NULL;
00051     m_computation_engine=NULL;
00052 
00053     SG_ADD((CSGObject**)&m_trace_sampler, "trace_sampler",
00054         "Trace sampler for the log operator", MS_NOT_AVAILABLE);
00055 
00056     SG_ADD((CSGObject**)&m_operator_log, "operator_log",
00057         "The log operator function", MS_NOT_AVAILABLE);
00058 
00059     SG_ADD((CSGObject**)&m_computation_engine, "computation_engine",
00060         "The computation engine for the jobs", MS_NOT_AVAILABLE);
00061 }
00062 
00063 CLogDetEstimator::~CLogDetEstimator()
00064 {
00065     SG_UNREF(m_trace_sampler);
00066     SG_UNREF(m_operator_log);
00067     SG_UNREF(m_computation_engine);
00068 }
00069 
00070 SGVector<float64_t> CLogDetEstimator::sample(index_t num_estimates)
00071 {
00072     SG_DEBUG("Entering\n");
00073     SG_INFO("Computing %d log-det estimates\n", num_estimates);
00074 
00075     REQUIRE(m_operator_log, "Operator function is NULL\n");
00076     // call the precompute of operator function to compute the prerequisites
00077     m_operator_log->precompute();
00078 
00079     REQUIRE(m_trace_sampler, "Trace sampler is NULL\n");
00080     // call the precompute of the sampler
00081     m_trace_sampler->precompute();
00082 
00083     REQUIRE(m_operator_log->get_operator()->get_dimension()\
00084         ==m_trace_sampler->get_dimension(),
00085         "Mismatch in dimensions of the operator and trace-sampler, %d vs %d!\n",
00086         m_operator_log->get_operator()->get_dimension(),
00087         m_trace_sampler->get_dimension());
00088 
00089     // for storing the aggregators that submit_jobs return
00090     CDynamicObjectArray* aggregators=new CDynamicObjectArray();
00091     index_t num_trace_samples=m_trace_sampler->get_num_samples();
00092 
00093     for (index_t i=0; i<num_estimates; ++i)
00094     {
00095         for (index_t j=0; j<num_trace_samples; ++j)
00096         {
00097             SG_INFO("Computing log-determinant trace sample %d/%d\n", j,
00098                     num_trace_samples);
00099 
00100             SG_DEBUG("Creating job for estimate %d, trace sample %d/%d\n", i, j,
00101                     num_trace_samples);
00102             // get the trace sampler vector
00103             SGVector<float64_t> s=m_trace_sampler->sample(j);
00104             // create jobs with the sample vector and store the aggregator
00105             CJobResultAggregator* agg=m_operator_log->submit_jobs(s);
00106             aggregators->append_element(agg);
00107             SG_UNREF(agg);
00108         }
00109     }
00110 
00111     REQUIRE(m_computation_engine, "Computation engine is NULL\n");
00112 
00113     // wait for all the jobs to be completed
00114     SG_INFO("Waiting for jobs to finish\n");
00115     m_computation_engine->wait_for_all();
00116     SG_INFO("All jobs finished, aggregating results\n");
00117 
00118     // the samples vector which stores the estimates with averaging
00119     SGVector<float64_t> samples(num_estimates);
00120     samples.zero();
00121 
00122     // use the aggregators to find the final result
00123     // use the same order as job submission to combine results
00124     int32_t num_aggregates=aggregators->get_num_elements();
00125     index_t idx_row=0;
00126     index_t idx_col=0;
00127     for (int32_t i=0; i<num_aggregates; ++i)
00128     {
00129         // this cast is safe due to above way of building the array
00130         CJobResultAggregator* agg=dynamic_cast<CJobResultAggregator*>
00131             (aggregators->get_element(i));
00132         ASSERT(agg);
00133 
00134         // call finalize on all the aggregators, cast is safe again
00135         agg->finalize();
00136         CScalarResult<float64_t>* r=dynamic_cast<CScalarResult<float64_t>*>
00137             (agg->get_final_result());
00138         ASSERT(r);
00139 
00140         // iterate through indices, group results in the same way as jobs
00141         samples[idx_col]+=r->get_result();
00142         idx_row++;
00143         if (idx_row>=num_trace_samples)
00144         {
00145             idx_row=0;
00146             idx_col++;
00147         }
00148 
00149         SG_UNREF(agg);
00150     }
00151 
00152     // clear all aggregators
00153     SG_UNREF(aggregators)
00154 
00155     SG_INFO("Finished computing %d log-det estimates\n", num_estimates);
00156 
00157     SG_DEBUG("Leaving\n");
00158     return samples;
00159 }
00160 
00161 SGMatrix<float64_t> CLogDetEstimator::sample_without_averaging(
00162     index_t num_estimates)
00163 {
00164     SG_DEBUG("Entering...\n")
00165 
00166     REQUIRE(m_operator_log, "Operator function is NULL\n");
00167     // call the precompute of operator function to compute all prerequisites
00168     m_operator_log->precompute();
00169 
00170     REQUIRE(m_trace_sampler, "Trace sampler is NULL\n");
00171     // call the precompute of the sampler
00172     m_trace_sampler->precompute();
00173 
00174     // for storing the aggregators that submit_jobs return
00175     CDynamicObjectArray aggregators;
00176     index_t num_trace_samples=m_trace_sampler->get_num_samples();
00177 
00178     for (index_t i=0; i<num_estimates; ++i)
00179     {
00180         for (index_t j=0; j<num_trace_samples; ++j)
00181         {
00182             // get the trace sampler vector
00183             SGVector<float64_t> s=m_trace_sampler->sample(j);
00184             // create jobs with the sample vector and store the aggregator
00185             CJobResultAggregator* agg=m_operator_log->submit_jobs(s);
00186             aggregators.append_element(agg);
00187             SG_UNREF(agg);
00188         }
00189     }
00190 
00191     REQUIRE(m_computation_engine, "Computation engine is NULL\n");
00192     // wait for all the jobs to be completed
00193     m_computation_engine->wait_for_all();
00194 
00195     // the samples matrix which stores the estimates without averaging
00196     // dimension: number of trace samples x number of log-det estimates
00197     SGMatrix<float64_t> samples(num_trace_samples, num_estimates);
00198 
00199     // use the aggregators to find the final result
00200     int32_t num_aggregates=aggregators.get_num_elements();
00201     for (int32_t i=0; i<num_aggregates; ++i)
00202     {
00203         CJobResultAggregator* agg=dynamic_cast<CJobResultAggregator*>
00204             (aggregators.get_element(i));
00205         if (!agg)
00206             SG_ERROR("Element is not CJobResultAggregator type!\n");
00207 
00208         // call finalize on all the aggregators
00209         agg->finalize();
00210         CScalarResult<float64_t>* r=dynamic_cast<CScalarResult<float64_t>*>
00211             (agg->get_final_result());
00212         if (!r)
00213             SG_ERROR("Result is not CScalarResult type!\n");
00214 
00215         // its important that we don't just unref the result here
00216         index_t idx_row=i%num_trace_samples;
00217         index_t idx_col=i/num_trace_samples;
00218         samples(idx_row, idx_col)=r->get_result();
00219         SG_UNREF(agg);
00220     }
00221 
00222     // clear all aggregators
00223     aggregators.clear_array();
00224 
00225     SG_DEBUG("Leaving\n")
00226     return samples;
00227 }
00228 
00229 }
00230 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

SHOGUN Machine Learning Toolbox - Documentation