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/config.h> 00011 00012 #ifdef HAVE_EIGEN3 00013 #include <shogun/base/Parameter.h> 00014 #include <shogun/lib/SGVector.h> 00015 #include <shogun/lib/SGMatrix.h> 00016 #include <shogun/mathematics/linalg/linsolver/LinearSolver.h> 00017 #include <shogun/mathematics/linalg/linop/DenseMatrixOperator.h> 00018 #include <shogun/mathematics/linalg/linop/SparseMatrixOperator.h> 00019 #include <shogun/mathematics/linalg/ratapprox/logdet/opfunc/LogRationalApproximationIndividual.h> 00020 #include <shogun/mathematics/linalg/ratapprox/logdet/computation/job/RationalApproximationIndividualJob.h> 00021 #include <shogun/mathematics/linalg/ratapprox/logdet/computation/aggregator/IndividualJobResultAggregator.h> 00022 #include <shogun/lib/computation/engine/IndependentComputationEngine.h> 00023 #include <typeinfo> 00024 00025 namespace shogun 00026 { 00027 00028 CLogRationalApproximationIndividual::CLogRationalApproximationIndividual() 00029 : CRationalApproximation(NULL, NULL, NULL, 0, OF_LOG) 00030 { 00031 init(); 00032 00033 SG_GCDEBUG("%s created (%p)\n", this->get_name(), this) 00034 } 00035 00036 CLogRationalApproximationIndividual::CLogRationalApproximationIndividual( 00037 CMatrixOperator<float64_t>* linear_operator, 00038 CIndependentComputationEngine* computation_engine, 00039 CEigenSolver* eigen_solver, 00040 CLinearSolver<complex128_t, float64_t>* linear_solver, 00041 float64_t desired_accuracy) 00042 : CRationalApproximation(linear_operator, computation_engine, 00043 eigen_solver, desired_accuracy, OF_LOG) 00044 { 00045 init(); 00046 00047 m_linear_solver=linear_solver; 00048 SG_REF(m_linear_solver); 00049 00050 SG_GCDEBUG("%s created (%p)\n", this->get_name(), this) 00051 } 00052 00053 void CLogRationalApproximationIndividual::init() 00054 { 00055 m_linear_solver=NULL; 00056 00057 SG_ADD((CSGObject**)&m_linear_solver, "linear_solver", 00058 "Linear solver for complex systems", MS_NOT_AVAILABLE); 00059 } 00060 00061 CLogRationalApproximationIndividual::~CLogRationalApproximationIndividual() 00062 { 00063 SG_UNREF(m_linear_solver); 00064 00065 SG_GCDEBUG("%s destroyed (%p)\n", this->get_name(), this) 00066 } 00067 00068 CJobResultAggregator* CLogRationalApproximationIndividual::submit_jobs( 00069 SGVector<float64_t> sample) 00070 { 00071 SG_DEBUG("OperatorFunction::submit_jobs(): Entering..\n"); 00072 REQUIRE(sample.vector, "Sample is not initialized!\n"); 00073 REQUIRE(m_linear_operator, "Operator is not initialized!\n"); 00074 REQUIRE(m_computation_engine, "Computation engine is NULL\n"); 00075 00076 // create the aggregator with sample, and the multiplier 00077 CIndividualJobResultAggregator* agg=new CIndividualJobResultAggregator( 00078 m_linear_operator, sample, m_constant_multiplier); 00079 // we don't want the aggregator to be destroyed when the job is unref-ed 00080 SG_REF(agg); 00081 00082 // this enum will save from repeated typechecking for all jobs 00083 enum typeID {DENSE=1, SPARSE, UNKNOWN} operator_type=UNKNOWN; 00084 00085 // create a complex copy of the matrix linear operator 00086 CMatrixOperator<complex128_t>* complex_op=NULL; 00087 if (typeid(*m_linear_operator)==typeid(CDenseMatrixOperator<float64_t>)) 00088 { 00089 operator_type=DENSE; 00090 00091 CDenseMatrixOperator<float64_t>* op 00092 =dynamic_cast<CDenseMatrixOperator<float64_t>*>(m_linear_operator); 00093 00094 REQUIRE(op->get_matrix_operator().matrix, "Matrix is not initialized!\n"); 00095 00096 // create complex dense matrix operator 00097 complex_op=static_cast<CDenseMatrixOperator<complex128_t>*>(*op); 00098 } 00099 else if (typeid(*m_linear_operator)==typeid(CSparseMatrixOperator<float64_t>)) 00100 { 00101 operator_type=SPARSE; 00102 00103 CSparseMatrixOperator<float64_t>* op 00104 =dynamic_cast<CSparseMatrixOperator<float64_t>*>(m_linear_operator); 00105 00106 REQUIRE(op->get_matrix_operator().sparse_matrix, "Matrix is not initialized!\n"); 00107 00108 // create complex sparse matrix operator 00109 complex_op=static_cast<CSparseMatrixOperator<complex128_t>*>(*op); 00110 } 00111 else 00112 { 00113 // something weird happened 00114 SG_ERROR("OperatorFunction::submit_jobs(): Unknown MatrixOperator given!\n"); 00115 } 00116 00117 // create num_shifts number of jobs for current sample vector 00118 for (index_t i=0; i<m_num_shifts; ++i) 00119 { 00120 // create a deep copy of the operator 00121 CMatrixOperator<complex128_t>* shifted_op=NULL; 00122 00123 switch(operator_type) 00124 { 00125 case DENSE: 00126 shifted_op=new CDenseMatrixOperator<complex128_t> 00127 (*dynamic_cast<CDenseMatrixOperator<complex128_t>*>(complex_op)); 00128 break; 00129 case SPARSE: 00130 shifted_op=new CSparseMatrixOperator<complex128_t> 00131 (*dynamic_cast<CSparseMatrixOperator<complex128_t>*>(complex_op)); 00132 break; 00133 default: 00134 break; 00135 } 00136 00137 REQUIRE(shifted_op, "OperatorFunction::submit_jobs():" 00138 "MatrixOperator typeinfo was not detected!\n"); 00139 00140 // move the shift inside the operator 00141 // (see CRationalApproximation) 00142 SGVector<complex128_t> diag=shifted_op->get_diagonal(); 00143 for (index_t j=0; j<diag.vlen; ++j) 00144 diag[j]-=m_shifts[i]; 00145 shifted_op->set_diagonal(diag); 00146 00147 // create a job and submit to the engine 00148 CRationalApproximationIndividualJob* job 00149 =new CRationalApproximationIndividualJob(agg, m_linear_solver, 00150 shifted_op, sample, m_weights[i]); 00151 SG_REF(job); 00152 00153 m_computation_engine->submit_job(job); 00154 00155 // we can safely unref the job here, computation engine takes it from here 00156 SG_UNREF(job); 00157 } 00158 00159 SG_UNREF(complex_op); 00160 00161 SG_DEBUG("OperatorFunction::submit_jobs(): Leaving..\n"); 00162 return agg; 00163 } 00164 00165 } 00166 #endif // HAVE_EIGEN3