SHOGUN
v3.2.0
|
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