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/lib/SGVector.h> 00014 #include <shogun/lib/computation/jobresult/ScalarResult.h> 00015 #include <shogun/mathematics/eigen3.h> 00016 #include <shogun/mathematics/linalg/linop/LinearOperator.h> 00017 #include <shogun/mathematics/linalg/linsolver/CGMShiftedFamilySolver.h> 00018 #include <shogun/mathematics/linalg/ratapprox/logdet/computation/job/RationalApproximationCGMJob.h> 00019 #include <shogun/base/Parameter.h> 00020 00021 using namespace Eigen; 00022 00023 namespace shogun 00024 { 00025 00026 CRationalApproximationCGMJob::CRationalApproximationCGMJob() 00027 : CIndependentJob() 00028 { 00029 init(); 00030 } 00031 00032 CRationalApproximationCGMJob::CRationalApproximationCGMJob( 00033 CStoreScalarAggregator<float64_t>* aggregator, 00034 CCGMShiftedFamilySolver* linear_solver, 00035 CLinearOperator<float64_t>* linear_operator, 00036 SGVector<float64_t> vector, 00037 SGVector<complex128_t> shifts, 00038 SGVector<complex128_t> weights, 00039 float64_t const_multiplier) 00040 : CIndependentJob((CJobResultAggregator*)aggregator) 00041 { 00042 init(); 00043 00044 m_linear_solver=linear_solver; 00045 SG_REF(m_linear_solver); 00046 00047 m_operator=linear_operator; 00048 SG_REF(m_operator); 00049 00050 m_vector=vector; 00051 00052 m_shifts=shifts; 00053 m_weights=weights; 00054 m_const_multiplier=const_multiplier; 00055 } 00056 00057 void CRationalApproximationCGMJob::init() 00058 { 00059 m_linear_solver=NULL; 00060 m_operator=NULL; 00061 m_const_multiplier=0.0; 00062 00063 SG_ADD((CSGObject**)&m_linear_solver, "linear_solver", 00064 "Linear solver for complex-shifted system", MS_NOT_AVAILABLE); 00065 00066 SG_ADD((CSGObject**)&m_operator, "linear_operator", 00067 "Linear operator", MS_NOT_AVAILABLE); 00068 00069 SG_ADD(&m_vector, "trace_sample", 00070 "Sample vector to apply linear operator on", MS_NOT_AVAILABLE); 00071 00072 SG_ADD(&m_weights, "complex_shifts", 00073 "Shifts in the linear systems to be solved", MS_NOT_AVAILABLE); 00074 00075 SG_ADD(&m_weights, "complex_weights", 00076 "Weights to be multiplied to the solution vector", MS_NOT_AVAILABLE); 00077 00078 SG_ADD(&m_const_multiplier, "constant_multiplier", 00079 "Constant multiplier to be multiplied with the final solution", MS_NOT_AVAILABLE); 00080 } 00081 00082 CRationalApproximationCGMJob::~CRationalApproximationCGMJob() 00083 { 00084 SG_UNREF(m_linear_solver); 00085 SG_UNREF(m_operator); 00086 } 00087 00088 void CRationalApproximationCGMJob::compute() 00089 { 00090 SG_DEBUG("Entering\n"); 00091 00092 REQUIRE(m_aggregator, "Job result aggregator is not set!\n"); 00093 REQUIRE(m_operator, "Operator is not set!\n"); 00094 REQUIRE(m_vector.vector, "Vector is not set!\n"); 00095 REQUIRE(m_shifts.vector, "Shifts are not set!\n"); 00096 REQUIRE(m_weights.vector, "Weights are not set!\n"); 00097 REQUIRE(m_operator->get_dimension()==m_vector.vlen, 00098 "Dimension mismatch! %d vs %d\n", m_operator->get_dimension(), m_vector.vlen); 00099 REQUIRE(m_shifts.vlen==m_weights.vlen, 00100 "Number of shifts and weights are not equal!\n"); 00101 00102 // solve the linear system with the sample vector 00103 SGVector<complex128_t> vec=m_linear_solver->solve_shifted_weighted( 00104 m_operator, m_vector, m_shifts, m_weights); 00105 00106 // Take negative (see CRationalApproximation for the formula) 00107 Map<VectorXcd> v(vec.vector, vec.vlen); 00108 v=-v; 00109 00110 // take out the imaginary part of the result before 00111 // applying linear operator 00112 SGVector<float64_t> agg=m_operator->apply(vec.get_imag()); 00113 00114 // perform dot product 00115 Map<VectorXd> map_agg(agg.vector, agg.vlen); 00116 Map<VectorXd> map_vector(m_vector.vector, m_vector.vlen); 00117 float64_t result=map_vector.dot(map_agg); 00118 00119 result*=m_const_multiplier; 00120 00121 // form the final result into a scalar result and submit to the aggregator 00122 CScalarResult<float64_t>* final_result=new CScalarResult<float64_t>(result); 00123 SG_REF(final_result); 00124 00125 m_aggregator->submit_result(final_result); 00126 00127 SG_UNREF(final_result); 00128 00129 SG_DEBUG("Leaving\n"); 00130 } 00131 00132 } 00133 #endif // HAVE_EIGEN3