Skip to content

Commit 8c5558f

Browse files
jcollfontjcollfont
committed
TilhonovSVD accepts the input of precomputed SVDs. It is not capable of using them yet
1 parent a9b56c6 commit 8c5558f

8 files changed

+303
-407
lines changed

‎src/Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithTikhonovSVD_impl.cc‎

Lines changed: 31 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -67,21 +67,37 @@ using namespace SCIRun::Core::Algorithms::Inverse;
6767
/////// prealocate Matrices for inverse compuation
6868
/// This function precalcualtes the SVD of the forward matrix and prepares singular vectors and values for posterior computations
6969
///////////////////////////////////////////////////////////////////
70-
void SolveInverseProblemWithTikhonovSVD_impl::preAlocateInverseMatrices(const SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_, const SCIRun::Core::Datatypes::DenseMatrix& measuredData_ , const SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_)
70+
void SolveInverseProblemWithTikhonovSVD_impl::preAlocateInverseMatrices(const SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_, const SCIRun::Core::Datatypes::DenseMatrix& measuredData_ , const SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& matrixU_, const SCIRun::Core::Datatypes::DenseMatrix& singularValues_, const SCIRun::Core::Datatypes::DenseMatrix& matrixV_)
7171
{
7272

73-
// Compute the SVD of the forward matrix
74-
Eigen::JacobiSVD<SCIRun::Core::Datatypes::DenseMatrix::EigenBase> svd( forwardMatrix_, Eigen::ComputeFullU | Eigen::ComputeFullV);
75-
SVDdecomposition = svd;
73+
// if (!matrixU_.size())&&(!matrixV_.size())&&(!singularValues_.size()){
74+
// std::cout << "Precomputed SVD variables found as an input" << std::endl;
75+
//
76+
// svd_MatrixU = *matrixU_;
77+
// svd_MatrixV = *matrixV_;
78+
// svd_SingularValues = *singularValues_;
79+
//
80+
// }else{
7681

77-
// determine rank
78-
rank = SVDdecomposition.nonzeroSingularValues();
82+
std::cout << "No precomputed SVD... computing now" << std::endl;
7983

80-
// Compute the projection of data y on the left singular vectors
81-
auto tempUy = SVDdecomposition.matrixU().transpose() * (measuredData_);
84+
// Compute the SVD of the forward matrix
85+
Eigen::JacobiSVD<SCIRun::Core::Datatypes::DenseMatrix::EigenBase> SVDdecomposition( forwardMatrix_, Eigen::ComputeFullU | Eigen::ComputeFullV);
8286

83-
Uy = tempUy;
87+
// alocate the left and right singular vectors and the singular values
88+
svd_MatrixU = SVDdecomposition.matrixU();
89+
svd_MatrixV = SVDdecomposition.matrixV();
90+
svd_SingularValues = SVDdecomposition.singularValues();
8491

92+
// determine rank
93+
rank = SVDdecomposition.nonzeroSingularValues();
94+
95+
// Compute the projection of data y on the left singular vectors
96+
auto tempUy = SVDdecomposition.matrixU().transpose() * (measuredData_);
97+
98+
Uy = tempUy;
99+
100+
// }
85101
}
86102

87103
//////////////////////////////////////////////////////////////////////
@@ -91,8 +107,8 @@ SCIRun::Core::Datatypes::DenseMatrix SolveInverseProblemWithTikhonovSVD_impl::co
91107
{
92108

93109
// prealocate matrices
94-
const int N = SVDdecomposition.matrixV().cols();
95-
const int M = SVDdecomposition.matrixU().rows();
110+
const int N = svd_MatrixV.cols();
111+
const int M = svd_MatrixU.rows();
96112
const int numTimeSamples = Uy.ncols();
97113
DenseMatrix solution(DenseMatrix::Zero(N,numTimeSamples));
98114
DenseMatrix tempInverse(DenseMatrix::Zero(N,M));
@@ -101,15 +117,15 @@ SCIRun::Core::Datatypes::DenseMatrix SolveInverseProblemWithTikhonovSVD_impl::co
101117
for (int rr=0; rr<rank ; rr++)
102118
{
103119
// evaluate filter factor
104-
double singVal = SVDdecomposition.singularValues()[rr];
120+
double singVal = svd_SingularValues[rr];
105121
double filterFactor_i = singVal / ( lambda_sq + singVal * singVal ) * Uy(rr);
106-
122+
107123
// u[date solution
108-
solution += filterFactor_i * SVDdecomposition.matrixV().col(rr);
124+
solution += filterFactor_i * svd_MatrixV.col(rr);
109125

110126
// update inverse operator
111127
if (inverseCalculation)
112-
tempInverse += filterFactor_i * ( SVDdecomposition.matrixV().col(rr) * SVDdecomposition.matrixU().col(rr).transpose() );
128+
tempInverse += filterFactor_i * ( svd_MatrixV.col(rr) * svd_MatrixU.col(rr).transpose() );
113129
}
114130

115131
// output solutions

‎src/Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithTikhonovSVD_impl.h‎

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -58,20 +58,23 @@ namespace Inverse {
5858

5959

6060
public:
61-
SolveInverseProblemWithTikhonovSVD_impl(const SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_, const SCIRun::Core::Datatypes::DenseMatrix& measuredData_ , const SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_)
61+
SolveInverseProblemWithTikhonovSVD_impl(const SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_, const SCIRun::Core::Datatypes::DenseMatrix& measuredData_ , const SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& matrixU_, const SCIRun::Core::Datatypes::DenseMatrix& singularValues_, const SCIRun::Core::Datatypes::DenseMatrix& matrixV_)
6262
{
63-
preAlocateInverseMatrices( forwardMatrix_, measuredData_ , sourceWeighting_, sensorWeighting_);
63+
preAlocateInverseMatrices( forwardMatrix_, measuredData_ , sourceWeighting_, sensorWeighting_, matrixU_, singularValues_, matrixV_ );
6464
};
6565

6666
private:
6767

6868
// Data Members
6969
int rank;
70-
Eigen::JacobiSVD<SCIRun::Core::Datatypes::DenseMatrix::EigenBase> SVDdecomposition;
70+
SCIRun::Core::Datatypes::DenseMatrix svd_MatrixU;
71+
SCIRun::Core::Datatypes::DenseColumnMatrix svd_SingularValues;
72+
SCIRun::Core::Datatypes::DenseMatrix svd_MatrixV;
73+
7174
SCIRun::Core::Datatypes::DenseMatrix Uy;
7275

7376
// Methods
74-
void preAlocateInverseMatrices(const SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_, const SCIRun::Core::Datatypes::DenseMatrix& measuredData_ , const SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_);
77+
void preAlocateInverseMatrices(const SCIRun::Core::Datatypes::DenseMatrix& forwardMatrix_, const SCIRun::Core::Datatypes::DenseMatrix& measuredData_ , const SCIRun::Core::Datatypes::DenseMatrix& sourceWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& sensorWeighting_, const SCIRun::Core::Datatypes::DenseMatrix& matrixU_, const SCIRun::Core::Datatypes::DenseMatrix& singularValues_, const SCIRun::Core::Datatypes::DenseMatrix& matrixV_);
7578

7679
virtual SCIRun::Core::Datatypes::DenseMatrix computeInverseSolution( double lambda_sq, bool inverseCalculation) const;
7780
// bool checkInputMatrixSizes(); // DEFINED IN PARENT, MIGHT WANT TO OVERRIDE SOME OTHER TIME

‎src/Core/Algorithms/Legacy/Inverse/TikhonovAlgoAbstractBase.cc‎

Lines changed: 13 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -60,11 +60,18 @@ using namespace SCIRun::Core::Logging;
6060
using namespace SCIRun::Core::Algorithms;
6161
using namespace SCIRun::Core::Algorithms::Inverse;
6262

63+
// shared inputs
6364
const AlgorithmInputName TikhonovAlgoAbstractBase::ForwardMatrix("ForwardMatrix");
6465
const AlgorithmInputName TikhonovAlgoAbstractBase::MeasuredPotentials("MeasuredPotentials");
6566
const AlgorithmInputName TikhonovAlgoAbstractBase::WeightingInSourceSpace("WeightingInSourceSpace");
6667
const AlgorithmInputName TikhonovAlgoAbstractBase::WeightingInSensorSpace("WeightingInSensorSpace");
6768

69+
// Inputs specific from the Tikhonov SVD module
70+
const AlgorithmInputName TikhonovAlgoAbstractBase::matrixU("matrixU");
71+
const AlgorithmInputName TikhonovAlgoAbstractBase::singularValues("singularValues");
72+
const AlgorithmInputName TikhonovAlgoAbstractBase::matrixV("matrixV");
73+
74+
// outputs
6875
const AlgorithmOutputName TikhonovAlgoAbstractBase::InverseSolution("InverseSolution");
6976
const AlgorithmOutputName TikhonovAlgoAbstractBase::RegularizationParameter("RegularizationParameter");
7077
const AlgorithmOutputName TikhonovAlgoAbstractBase::RegInverse("RegInverse");
@@ -224,7 +231,12 @@ AlgorithmOutput TikhonovAlgoAbstractBase::run(const AlgorithmInput & input) cons
224231
int regularizationSolutionSubcase_ = get(Parameters::regularizationSolutionSubcase).toInt();
225232
int regularizationResidualSubcase_ = get(Parameters::regularizationResidualSubcase).toInt();
226233

227-
algoImpl = new SolveInverseProblemWithTikhonovSVD_impl(*castMatrix::toDense(forwardMatrix_), *castMatrix::toDense(measuredData_), *castMatrix::toDense(sourceWeighting_), *castMatrix::toDense(sensorWeighting_));
234+
// get TikhonovSVD special inputs
235+
auto matrixU_ = input.get<Matrix>(TikhonovAlgoAbstractBase::matrixU);
236+
auto singularValues_ = input.get<Matrix>(TikhonovAlgoAbstractBase::singularValues);
237+
auto matrixV_ = input.get<Matrix>(TikhonovAlgoAbstractBase::matrixV);
238+
239+
algoImpl = new SolveInverseProblemWithTikhonovSVD_impl( *castMatrix::toDense(forwardMatrix_), *castMatrix::toDense(measuredData_), *castMatrix::toDense(sourceWeighting_), *castMatrix::toDense(sensorWeighting_), *castMatrix::toDense(matrixU_), *castMatrix::toDense(singularValues_), *castMatrix::toDense(matrixV_) );
228240
}
229241
else if ( TikhonovImplementation_gotten == std::string("TikhonovTSVD") ){
230242
THROW_ALGORITHM_PROCESSING_ERROR("Tikhonov TSVD not implemented yet");

‎src/Core/Algorithms/Legacy/Inverse/TikhonovAlgoAbstractBase.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -65,6 +65,11 @@ namespace Inverse {
6565
static const AlgorithmInputName WeightingInSourceSpace;
6666
static const AlgorithmInputName WeightingInSensorSpace;
6767

68+
// input names particular for the SVD case
69+
static const AlgorithmInputName matrixU;
70+
static const AlgorithmInputName singularValues;
71+
static const AlgorithmInputName matrixV;
72+
6873
// define output names
6974
static const AlgorithmOutputName InverseSolution;
7075
static const AlgorithmOutputName RegularizationParameter;

‎src/Interface/Modules/Inverse/SolveInverseProblemWithTikhonovSVDDialog.cc‎

Lines changed: 0 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -49,8 +49,6 @@ SolveInverseProblemWithTikhonovSVDDialog::SolveInverseProblemWithTikhonovSVDDial
4949
lambdaMethod_.insert(StringPair("Slider", "slider"));
5050
lambdaMethod_.insert(StringPair("L-curve", "lcurve"));
5151

52-
WidgetStyleMixin::tabStyle(inputTabWidget_);
53-
5452
addDoubleLineEditManager(lCurveLambdaLineEdit_, Parameters::LambdaCorner);
5553
addSpinBoxManager(lambdaNumberSpinBox_, Parameters::LambdaNum);
5654
addDoubleSpinBoxManager(lambdaDoubleSpinBox_, Parameters::LambdaFromDirectEntry);
@@ -62,10 +60,6 @@ SolveInverseProblemWithTikhonovSVDDialog::SolveInverseProblemWithTikhonovSVDDial
6260

6361
addDoubleSpinBoxManager(lambdaSliderDoubleSpinBox_, Parameters::LambdaSliderValue);
6462

65-
addRadioButtonGroupManager({ autoRadioButton_, underRadioButton_, overRadioButton_ }, Parameters::regularizationChoice);
66-
addRadioButtonGroupManager({ solutionConstraintRadioButton_, squaredSolutionRadioButton_ }, Parameters::regularizationSolutionSubcase);
67-
addRadioButtonGroupManager({ residualConstraintRadioButton_, squaredResidualSolutionRadioButton_ }, Parameters::regularizationResidualSubcase);
68-
6963
addComboBoxManager(lambdaMethodComboBox_, Parameters::RegularizationMethod, lambdaMethod_);
7064
addTextEditManager(lCurveTextEdit_, Parameters::LCurveText);
7165

0 commit comments

Comments
 (0)