Skip to content

Commit

Permalink
Refactoring (#4103)
Browse files Browse the repository at this point in the history
* Refactoring

* Update LogDetEstimator.cpp

* rebase and edits

* fixes RationalApproximation unit test

* style fixes
  • Loading branch information
shubham808 authored and karlnapf committed Apr 2, 2018
1 parent bccd76c commit 8d578f4
Show file tree
Hide file tree
Showing 21 changed files with 193 additions and 576 deletions.
7 changes: 2 additions & 5 deletions examples/undocumented/python/mathematics_logdet.py
Expand Up @@ -27,7 +27,6 @@ def mathematics_logdet (matrix=mtx,max_iter_eig=1000,max_iter_lin=1000,num_sampl
from shogun.Mathematics import ProbingSampler
from shogun.Mathematics import LogDetEstimator
from shogun.Mathematics import Statistics
from shogun.Library import SerialComputationEngine

# creating the linear operator, eigen-solver
op=RealSparseMatrixOperator(A.tocsc())
Expand All @@ -50,16 +49,14 @@ def mathematics_logdet (matrix=mtx,max_iter_eig=1000,max_iter_lin=1000,num_sampl
lin_solver=CGMShiftedFamilySolver()
lin_solver.set_iteration_limit(max_iter_lin)

# computation engine
engine=SerialComputationEngine()

# set the desired accuracy tighter to obtain better results
# this determines the number of contour points in conformal mapping of
# the rational approximation of the Cauchy's integral of f(A)*s, f=log
desired_accuracy=1E-5

# creating the log-linear-operator function
op_func=LogRationalApproximationCGM(op, engine, eig_solver, lin_solver,\
op_func=LogRationalApproximationCGM(op, eig_solver, lin_solver,\
desired_accuracy)

# set the trace sampler to be probing sampler, in which samples are obtained
Expand All @@ -68,7 +65,7 @@ def mathematics_logdet (matrix=mtx,max_iter_eig=1000,max_iter_lin=1000,num_sampl
trace_sampler=ProbingSampler(op)

# estimating log-det
log_det_estimator=LogDetEstimator(trace_sampler, op_func, engine)
log_det_estimator=LogDetEstimator(trace_sampler, op_func)

# set the number of samples as required
estimates=log_det_estimator.sample(num_samples)
Expand Down
1 change: 0 additions & 1 deletion src/shogun/base/class_list.h
Expand Up @@ -69,7 +69,6 @@ namespace shogun {
*
*/
std::set<std::string> available_objects();

}

#endif /* __SG_CLASS_LIST_H__ */
7 changes: 3 additions & 4 deletions src/shogun/labels/BinaryLabels.cpp
Expand Up @@ -151,8 +151,7 @@ CLabels* CBinaryLabels::shallow_subset_copy()
return shallow_copy_labels;
}

CBinaryLabels::CBinaryLabels(const CDenseLabels& dense):
CDenseLabels(dense)
CBinaryLabels::CBinaryLabels(const CDenseLabels& dense) : CDenseLabels(dense)
{
ensure_valid();
}
Expand All @@ -179,8 +178,8 @@ namespace shogun
catch (const ShogunException& e)
{
SG_SERROR(
"Cannot convert %s to binary labels: %s\n",
orig->get_name(), e.what());
"Cannot convert %s to binary labels: %s\n", orig->get_name(),
e.what());
}

return Some<CBinaryLabels>::from_raw(nullptr);
Expand Down
6 changes: 2 additions & 4 deletions src/shogun/labels/DenseLabels.cpp
Expand Up @@ -32,14 +32,12 @@ CDenseLabels::CDenseLabels(int32_t num_lab)
linalg::zero(m_current_values);
}

CDenseLabels::CDenseLabels(const CDenseLabels& orig):
CLabels(orig),
m_labels(orig.m_labels)
CDenseLabels::CDenseLabels(const CDenseLabels& orig)
: CLabels(orig), m_labels(orig.m_labels)
{
init();
}


CDenseLabels::CDenseLabels(CFile* loader)
: CLabels()
{
Expand Down
3 changes: 1 addition & 2 deletions src/shogun/labels/Labels.cpp
Expand Up @@ -21,8 +21,7 @@ CLabels::CLabels()
}

CLabels::CLabels(const CLabels& orig)
: CSGObject(orig),
m_current_values(orig.m_current_values)
: CSGObject(orig), m_current_values(orig.m_current_values)
{
init();
}
Expand Down
3 changes: 2 additions & 1 deletion src/shogun/labels/MulticlassLabels.cpp
Expand Up @@ -215,7 +215,8 @@ namespace shogun
switch (orig->get_label_type())
{
case LT_MULTICLASS:
return Some<CMulticlassLabels>::from_raw(orig->as<CMulticlassLabels>());
return Some<CMulticlassLabels>::from_raw(
orig->as<CMulticlassLabels>());
case LT_DENSE_GENERIC:
return to_multiclass(orig->as<CDenseLabels>());
default:
Expand Down
1 change: 0 additions & 1 deletion src/shogun/machine/Machine.h
Expand Up @@ -455,6 +455,5 @@ class CMachine : public CSGObject
/** Mutex used to pause threads */
std::mutex m_mutex;
};

}
#endif // _MACHINE_H__
182 changes: 41 additions & 141 deletions src/shogun/mathematics/linalg/ratapprox/logdet/LogDetEstimator.cpp
Expand Up @@ -3,26 +3,21 @@
*
* Authors: Sunil Mahendrakar, Heiko Strathmann, Soumyajit De, Björn Esser
*/

#include <shogun/lib/common.h>
#include <shogun/lib/SGVector.h>
#include <shogun/base/Parallel.h>
#include <shogun/lib/SGMatrix.h>
#include <shogun/lib/DynamicObjectArray.h>
#include <shogun/lib/computation/engine/IndependentComputationEngine.h>
#include <shogun/lib/computation/engine/SerialComputationEngine.h>
#include <shogun/lib/computation/jobresult/ScalarResult.h>
#include <shogun/lib/computation/aggregator/JobResultAggregator.h>
#include <shogun/lib/SGVector.h>
#include <shogun/lib/common.h>
#include <shogun/mathematics/linalg/eigsolver/LanczosEigenSolver.h>
#include <shogun/mathematics/linalg/linop/DenseMatrixOperator.h>
#include <shogun/mathematics/linalg/linop/SparseMatrixOperator.h>
#include <shogun/mathematics/linalg/eigsolver/LanczosEigenSolver.h>
#include <shogun/mathematics/linalg/linsolver/CGMShiftedFamilySolver.h>
#include <shogun/mathematics/linalg/ratapprox/tracesampler/TraceSampler.h>
#include <shogun/mathematics/linalg/ratapprox/tracesampler/ProbingSampler.h>
#include <shogun/mathematics/linalg/ratapprox/tracesampler/NormalSampler.h>
#include <shogun/mathematics/linalg/ratapprox/opfunc/OperatorFunction.h>
#include <shogun/mathematics/linalg/ratapprox/logdet/LogDetEstimator.h>
#include <shogun/mathematics/linalg/ratapprox/logdet/opfunc/DenseMatrixExactLog.h>
#include <shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationCGM.h>
#include <shogun/mathematics/linalg/ratapprox/opfunc/OperatorFunction.h>
#include <shogun/mathematics/linalg/ratapprox/tracesampler/NormalSampler.h>
#include <shogun/mathematics/linalg/ratapprox/tracesampler/ProbingSampler.h>
#include <shogun/mathematics/linalg/ratapprox/tracesampler/TraceSampler.h>

namespace shogun
{
Expand All @@ -39,9 +34,6 @@ CLogDetEstimator::CLogDetEstimator(SGSparseMatrix<float64_t> sparse_mat)
{
init();

m_computation_engine=new CSerialComputationEngine();
SG_REF(m_computation_engine);

CSparseMatrixOperator<float64_t>* op=
new CSparseMatrixOperator<float64_t>(sparse_mat);

Expand All @@ -50,8 +42,8 @@ CLogDetEstimator::CLogDetEstimator(SGSparseMatrix<float64_t> sparse_mat)
CLanczosEigenSolver* eig_solver=new CLanczosEigenSolver(op);
CCGMShiftedFamilySolver* linear_solver=new CCGMShiftedFamilySolver();

m_operator_log=new CLogRationalApproximationCGM(op,m_computation_engine,
eig_solver,linear_solver,accuracy);
m_operator_log = new CLogRationalApproximationCGM(
op, eig_solver, linear_solver, accuracy);
SG_REF(m_operator_log);

#ifdef HAVE_COLPACK
Expand All @@ -62,15 +54,14 @@ CLogDetEstimator::CLogDetEstimator(SGSparseMatrix<float64_t> sparse_mat)

SG_REF(m_trace_sampler);

SG_INFO("LogDetEstimator:Using %s, %s with 1E-5 accuracy, %s as default\n",
m_computation_engine->get_name(), m_operator_log->get_name(),
m_trace_sampler->get_name());
SG_INFO(
"LogDetEstimator: %s with 1E-5 accuracy, %s as default\n",
m_operator_log->get_name(), m_trace_sampler->get_name());
}
#endif //HAVE_LAPACK

CLogDetEstimator::CLogDetEstimator(CTraceSampler* trace_sampler,
COperatorFunction<float64_t>* operator_log,
CIndependentComputationEngine* computation_engine)
CLogDetEstimator::CLogDetEstimator(
CTraceSampler* trace_sampler, COperatorFunction<float64_t>* operator_log)
: CSGObject()
{
init();
Expand All @@ -80,32 +71,24 @@ CLogDetEstimator::CLogDetEstimator(CTraceSampler* trace_sampler,

m_operator_log=operator_log;
SG_REF(m_operator_log);

m_computation_engine=computation_engine;
SG_REF(m_computation_engine);
}

void CLogDetEstimator::init()
{
m_trace_sampler=NULL;
m_operator_log=NULL;
m_computation_engine=NULL;

SG_ADD((CSGObject**)&m_trace_sampler, "trace_sampler",
"Trace sampler for the log operator", MS_NOT_AVAILABLE);

SG_ADD((CSGObject**)&m_operator_log, "operator_log",
"The log operator function", MS_NOT_AVAILABLE);

SG_ADD((CSGObject**)&m_computation_engine, "computation_engine",
"The computation engine for the jobs", MS_NOT_AVAILABLE);
}

CLogDetEstimator::~CLogDetEstimator()
{
SG_UNREF(m_trace_sampler);
SG_UNREF(m_operator_log);
SG_UNREF(m_computation_engine);
}

CTraceSampler* CLogDetEstimator::get_trace_sampler(void) const
Expand All @@ -114,12 +97,6 @@ CTraceSampler* CLogDetEstimator::get_trace_sampler(void) const
return m_trace_sampler;
}

CIndependentComputationEngine* CLogDetEstimator::get_computation_engine(void) const
{
SG_REF(m_computation_engine);
return m_computation_engine;
}

COperatorFunction<float64_t>* CLogDetEstimator::get_operator_function(void) const
{
SG_REF(m_operator_log);
Expand All @@ -145,71 +122,26 @@ SGVector<float64_t> CLogDetEstimator::sample(index_t num_estimates)
m_operator_log->get_operator()->get_dimension(),
m_trace_sampler->get_dimension());

// for storing the aggregators that submit_jobs return
CDynamicObjectArray* aggregators=new CDynamicObjectArray();
index_t num_trace_samples=m_trace_sampler->get_num_samples();

for (index_t i=0; i<num_estimates; ++i)
{
for (index_t j=0; j<num_trace_samples; ++j)
{
SG_INFO("Computing log-determinant trace sample %d/%d\n", j,
num_trace_samples);

SG_DEBUG("Creating job for estimate %d, trace sample %d/%d\n", i, j,
num_trace_samples);
// get the trace sampler vector
SGVector<float64_t> s=m_trace_sampler->sample(j);
// create jobs with the sample vector and store the aggregator
CJobResultAggregator* agg=m_operator_log->submit_jobs(s);
aggregators->append_element(agg);
SG_UNREF(agg);
}
}

REQUIRE(m_computation_engine, "Computation engine is NULL\n");

// wait for all the jobs to be completed
SG_INFO("Waiting for jobs to finish\n");
m_computation_engine->wait_for_all();
SG_INFO("All jobs finished, aggregating results\n");

// the samples vector which stores the estimates with averaging
SGVector<float64_t> samples(num_estimates);
samples.zero();

// use the aggregators to find the final result
// use the same order as job submission to combine results
int32_t num_aggregates=aggregators->get_num_elements();
index_t idx_row=0;
index_t idx_col=0;
for (int32_t i=0; i<num_aggregates; ++i)
for (index_t i = 0; i < num_estimates; ++i)
{
// this cast is safe due to above way of building the array
CJobResultAggregator* agg=dynamic_cast<CJobResultAggregator*>
(aggregators->get_element(i));
ASSERT(agg);

// call finalize on all the aggregators, cast is safe again
agg->finalize();
CScalarResult<float64_t>* r=dynamic_cast<CScalarResult<float64_t>*>
(agg->get_final_result());
ASSERT(r);

// iterate through indices, group results in the same way as jobs
samples[idx_col]+=r->get_result();
idx_row++;
if (idx_row>=num_trace_samples)
for (index_t j = 0; j < num_trace_samples; ++j)
{
idx_row=0;
idx_col++;
SG_INFO(
"Computing log-determinant trace sample %d/%d\n", j,
num_trace_samples);
// get the trace sampler vector
SGVector<float64_t> s = m_trace_sampler->sample(j);
// calculate the result for sample s
float64_t result = m_operator_log->compute(s);
{
samples[i] += result;
}
}
}

SG_UNREF(agg);
}

// clear all aggregators
SG_UNREF(aggregators)

SG_INFO("Finished computing %d log-det estimates\n", num_estimates);

Expand All @@ -230,60 +162,28 @@ SGMatrix<float64_t> CLogDetEstimator::sample_without_averaging(
// call the precompute of the sampler
m_trace_sampler->precompute();

// for storing the aggregators that submit_jobs return
CDynamicObjectArray aggregators;
index_t num_trace_samples=m_trace_sampler->get_num_samples();
index_t num_trace_samples = m_trace_sampler->get_num_samples();
SGMatrix<float64_t> samples(num_trace_samples, num_estimates);

for (index_t i=0; i<num_estimates; ++i)
for (index_t i = 0; i < num_estimates; ++i)
{
for (index_t j=0; j<num_trace_samples; ++j)
for (index_t j = 0; j < num_trace_samples; ++j)
{
SG_INFO(
"Computing log-determinant trace sample %d/%d\n", j,
num_trace_samples);
// get the trace sampler vector
SGVector<float64_t> s=m_trace_sampler->sample(j);
// create jobs with the sample vector and store the aggregator
CJobResultAggregator* agg=m_operator_log->submit_jobs(s);
aggregators.append_element(agg);
SG_UNREF(agg);
SGVector<float64_t> s = m_trace_sampler->sample(j);
// solve the result for s
float64_t result = m_operator_log->compute(s);
{
samples(i, j) = result;
}
}
}

REQUIRE(m_computation_engine, "Computation engine is NULL\n");
// wait for all the jobs to be completed
m_computation_engine->wait_for_all();

// the samples matrix which stores the estimates without averaging
// dimension: number of trace samples x number of log-det estimates
SGMatrix<float64_t> samples(num_trace_samples, num_estimates);

// use the aggregators to find the final result
int32_t num_aggregates=aggregators.get_num_elements();
for (int32_t i=0; i<num_aggregates; ++i)
{
CJobResultAggregator* agg=dynamic_cast<CJobResultAggregator*>
(aggregators.get_element(i));
if (!agg)
SG_ERROR("Element is not CJobResultAggregator type!\n");

// call finalize on all the aggregators
agg->finalize();
CScalarResult<float64_t>* r=dynamic_cast<CScalarResult<float64_t>*>
(agg->get_final_result());
if (!r)
SG_ERROR("Result is not CScalarResult type!\n");

// its important that we don't just unref the result here
index_t idx_row=i%num_trace_samples;
index_t idx_col=i/num_trace_samples;
samples(idx_row, idx_col)=r->get_result();
SG_UNREF(agg);
}

// clear all aggregators
aggregators.clear_array();

SG_DEBUG("Leaving\n")
return samples;
}

}

0 comments on commit 8d578f4

Please sign in to comment.