Skip to content

Commit 4ec0481

Browse files
jcollfontjcollfont
committed
TSVD implemented. The L-curve is not doing what I would like it to do. The scale in which it should select the lambdas must be linear and now it is linear on the exponent…. not sure how to fix this or how much it matters
1 parent 2507820 commit 4ec0481

8 files changed

+60
-37
lines changed

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

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -57,7 +57,7 @@ using namespace SCIRun::Core::Algorithms::Inverse;
5757

5858
/////////////////////////
5959
///////// compute Inverse solution
60-
DenseMatrix SolveInverseProblemWithStandardTikhonovImpl::computeInverseSolution( double lambda_sq, bool inverseCalculation) const
60+
DenseMatrix SolveInverseProblemWithStandardTikhonovImpl::computeInverseSolution( double lambda, bool inverseCalculation) const
6161
{
6262
//............................
6363
// OPERATIONS PERFORMED IN THIS SECTION:
@@ -79,7 +79,7 @@ using namespace SCIRun::Core::Algorithms::Inverse;
7979
DenseMatrix solution(sizeSolution,numTimeSamples);
8080
DenseMatrix G;
8181

82-
G = M1 + lambda_sq * M2;
82+
G = M1 + lambda * lambda * M2;
8383

8484
b = G.lu().solve(y).eval();
8585

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -68,7 +68,7 @@ namespace SCIRun {
6868

6969
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 int regularizationChoice_, const int regularizationSolutionSubcase_, const int regularizationResidualSubcase_ );
7070

71-
virtual SCIRun::Core::Datatypes::DenseMatrix computeInverseSolution( double lambda_sq, bool inverseCalculation) const;
71+
virtual SCIRun::Core::Datatypes::DenseMatrix computeInverseSolution( double lambda, bool inverseCalculation) const;
7272
// bool checkInputMatrixSizes(); // DEFINED IN PARENT, MIGHT WANT TO OVERRIDE SOME OTHER TIME
7373

7474
};

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

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -110,7 +110,7 @@ void SolveInverseProblemWithTikhonovSVD_impl::preAlocateInverseMatrices(const SC
110110
//////////////////////////////////////////////////////////////////////
111111
// THIS FUNCTION returns regularized solution by tikhonov method
112112
//////////////////////////////////////////////////////////////////////
113-
SCIRun::Core::Datatypes::DenseMatrix SolveInverseProblemWithTikhonovSVD_impl::computeInverseSolution( double lambda_sq, bool inverseCalculation ) const
113+
SCIRun::Core::Datatypes::DenseMatrix SolveInverseProblemWithTikhonovSVD_impl::computeInverseSolution( double lambda, bool inverseCalculation ) const
114114
{
115115

116116
// prealocate matrices
@@ -125,7 +125,7 @@ SCIRun::Core::Datatypes::DenseMatrix SolveInverseProblemWithTikhonovSVD_impl::co
125125
{
126126
// evaluate filter factor
127127
double singVal = svd_SingularValues[rr];
128-
double filterFactor_i = singVal / ( lambda_sq + singVal * singVal ) * Uy(rr);
128+
double filterFactor_i = singVal / ( lambda * lambda + singVal * singVal ) * Uy(rr);
129129

130130
// u[date solution
131131
solution += filterFactor_i * svd_MatrixV.col(rr);

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -82,7 +82,7 @@ namespace Inverse {
8282
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_);
8383
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_);
8484

85-
virtual SCIRun::Core::Datatypes::DenseMatrix computeInverseSolution( double lambda_sq, bool inverseCalculation) const;
85+
virtual SCIRun::Core::Datatypes::DenseMatrix computeInverseSolution( double lambda, bool inverseCalculation) const;
8686
// bool checkInputMatrixSizes(); // DEFINED IN PARENT, MIGHT WANT TO OVERRIDE SOME OTHER TIME
8787

8888

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

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -87,6 +87,7 @@ void SolveInverseProblemWithTikhonovTSVD_impl::preAlocateInverseMatrices(const S
8787

8888
// determine rank
8989
rank = svd_SingularValues.nrows();
90+
9091
}
9192

9293
void SolveInverseProblemWithTikhonovTSVD_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_)
@@ -110,7 +111,7 @@ void SolveInverseProblemWithTikhonovTSVD_impl::preAlocateInverseMatrices(const S
110111
//////////////////////////////////////////////////////////////////////
111112
// THIS FUNCTION returns regularized solution by tikhonov method
112113
//////////////////////////////////////////////////////////////////////
113-
SCIRun::Core::Datatypes::DenseMatrix SolveInverseProblemWithTikhonovTSVD_impl::computeInverseSolution( double truncationPoint, bool inverseCalculation ) const
114+
SCIRun::Core::Datatypes::DenseMatrix SolveInverseProblemWithTikhonovTSVD_impl::computeInverseSolution( double lambda, bool inverseCalculation ) const
114115
{
115116

116117
// prealocate matrices
@@ -120,6 +121,8 @@ SCIRun::Core::Datatypes::DenseMatrix SolveInverseProblemWithTikhonovTSVD_impl::c
120121
DenseMatrix solution(DenseMatrix::Zero(N,numTimeSamples));
121122
DenseMatrix tempInverse(DenseMatrix::Zero(N,M));
122123

124+
const int truncationPoint = Min( int(lambda), rank, int(9999999999999) );
125+
123126
// Compute inverse SolveInverseProblemWithTikhonovTSVD
124127
for (int rr=0; rr < truncationPoint ; rr++)
125128
{

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

Lines changed: 19 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -91,6 +91,7 @@ ALGORITHM_PARAMETER_DEF( Inverse, LCurveText);
9191
ALGORITHM_PARAMETER_DEF( Inverse, regularizationSolutionSubcase);
9292
ALGORITHM_PARAMETER_DEF( Inverse, regularizationResidualSubcase);
9393

94+
9495
TikhonovAlgoAbstractBase::TikhonovAlgoAbstractBase()
9596
{
9697
using namespace Parameters;
@@ -208,7 +209,7 @@ AlgorithmOutput TikhonovAlgoAbstractBase::run(const AlgorithmInput & input) cons
208209

209210
// Alocate Variable
210211
DenseMatrix solution;
211-
double lambda_sq = 0;
212+
double lambda = 0;
212213
double lambda_ = 0;
213214

214215
// check input MATRICES
@@ -224,6 +225,9 @@ AlgorithmOutput TikhonovAlgoAbstractBase::run(const AlgorithmInput & input) cons
224225
int regularizationResidualSubcase_ = get(Parameters::regularizationResidualSubcase).toInt();
225226

226227
algoImpl = new SolveInverseProblemWithStandardTikhonovImpl( *castMatrix::toDense(forwardMatrix_), *castMatrix::toDense(measuredData_), *castMatrix::toDense(sourceWeighting_), *castMatrix::toDense(sensorWeighting_), regularizationChoice_, regularizationSolutionSubcase_, regularizationResidualSubcase_);
228+
229+
230+
227231
}
228232
else if ( TikhonovImplementation_gotten == std::string("TikhonovSVD") ){
229233

@@ -237,12 +241,15 @@ AlgorithmOutput TikhonovAlgoAbstractBase::run(const AlgorithmInput & input) cons
237241
auto singularValues_ = input.get<Matrix>(TikhonovAlgoAbstractBase::singularValues);
238242
auto matrixV_ = input.get<Matrix>(TikhonovAlgoAbstractBase::matrixV);
239243

244+
240245
// If there is a missing matrix from the precomputed SVD input
241246
if ( (matrixU_ == NULL) || (singularValues_ == NULL) || ( matrixV_ == NULL) )
242247
algoImpl = new SolveInverseProblemWithTikhonovSVD_impl( *castMatrix::toDense(forwardMatrix_), *castMatrix::toDense(measuredData_), *castMatrix::toDense(sourceWeighting_), *castMatrix::toDense(sensorWeighting_) );
243248
else
244249
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_) );
245250

251+
252+
246253
}
247254
else if ( TikhonovImplementation_gotten == std::string("TikhonovTSVD") ){
248255

@@ -256,17 +263,23 @@ AlgorithmOutput TikhonovAlgoAbstractBase::run(const AlgorithmInput & input) cons
256263
auto singularValues_ = input.get<Matrix>(TikhonovAlgoAbstractBase::singularValues);
257264
auto matrixV_ = input.get<Matrix>(TikhonovAlgoAbstractBase::matrixV);
258265

266+
259267
// If there is a missing matrix from the precomputed SVD input
260268
if ( (matrixU_ == NULL) || (singularValues_ == NULL) || ( matrixV_ == NULL) )
261269
algoImpl = new SolveInverseProblemWithTikhonovTSVD_impl( *castMatrix::toDense(forwardMatrix_), *castMatrix::toDense(measuredData_), *castMatrix::toDense(sourceWeighting_), *castMatrix::toDense(sensorWeighting_) );
262270
else
263271
algoImpl = new SolveInverseProblemWithTikhonovTSVD_impl( *castMatrix::toDense(forwardMatrix_), *castMatrix::toDense(measuredData_), *castMatrix::toDense(sourceWeighting_), *castMatrix::toDense(sensorWeighting_), *castMatrix::toDense(matrixU_), *castMatrix::toDense(singularValues_), *castMatrix::toDense(matrixV_) );
272+
273+
274+
264275
}
265276
else{
266277
THROW_ALGORITHM_PROCESSING_ERROR("Not a valid Tikhonov Implementation selection");
278+
267279
}
268280

269-
// preAlocateInverseMatrices(forwardMatrix_,measuredData_,sourceWeighting_,sensorWeighting_);
281+
282+
270283

271284
//Get Regularization parameter(s) : Lambda
272285
if ((RegularizationMethod_gotten == "single") || (RegularizationMethod_gotten == "slider"))
@@ -292,21 +305,12 @@ AlgorithmOutput TikhonovAlgoAbstractBase::run(const AlgorithmInput & input) cons
292305
}
293306

294307
std::cout << "Lambda: " << lambda_ << std::endl;
295-
lambda_sq = lambda_ * lambda_;
308+
lambda = lambda_;
296309

297310

298-
// compute inverse solution
299-
solution = algoImpl->computeInverseSolution(lambda_sq, true);
311+
// compute final inverse solution
312+
solution = algoImpl->computeInverseSolution(lambda, true);
300313

301-
//
302-
// // set final result
303-
// inverseSolution_.reset(new DenseMatrix(solution));
304-
305-
// output regularization parameter
306-
// DenseColumnMatrix tempLambda(1);
307-
// tempLambda[0] = lambda_;
308-
//
309-
// regularizationParameter_. reset( new DenseColumnMatrix(tempLambda) );
310314

311315

312316
// Set outputs
@@ -337,7 +341,6 @@ double TikhonovAlgoAbstractBase::computeLcurve( const SCIRun::Core::Algorithms::
337341
const double lam_step = (log10(lambdaMax_) - log10(lambdaMin_)) / (nLambda-1);
338342
std::cout << "Lambda power step: " << lam_step << ". Number: "<< nLambda <<". Lambda min: " << lambdaMin_ << ". Lambda max: "<< lambdaMax_<< ". Ratio: "<< lambdaMax_ / lambdaMin_ << std::endl;
339343
double lambda = 0;
340-
double lambda_sq;
341344

342345
// prealocate vector of lambdas and eta and rho
343346
std::vector<double> lambdaArray(nLambda, 0.0);
@@ -357,11 +360,8 @@ double TikhonovAlgoAbstractBase::computeLcurve( const SCIRun::Core::Algorithms::
357360
lambdaArray[j] = lambdaArray[j-1] * pow(10.0,lam_step);
358361
}
359362

360-
// set current lambda
361-
lambda_sq = lambdaArray[j] * lambdaArray[j];
362-
363363
// COMPUTE INVERSE SOLUTION
364-
solution = algoImpl->computeInverseSolution( lambda_sq, false);
364+
solution = algoImpl->computeInverseSolution( lambdaArray[j], false);
365365

366366

367367
// if using source regularization matrix, apply it to compute Rx (for the eta computations)

‎src/Modules/Legacy/Inverse/SolveInverseProblemWithTikhonovSVD.cc‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -115,7 +115,6 @@ void SolveInverseProblemWithTikhonovSVD::execute()
115115
setAlgoIntFromState(Parameters::LambdaCorner);
116116
setAlgoStringFromState(Parameters::LCurveText);
117117

118-
// run
119118
// run
120119
auto output = algo().run( withInputData((ForwardMatrix, forward_matrix_h)(MeasuredPotentials,hMatrixMeasDat)(MeasuredPotentials,hMatrixMeasDat)(WeightingInSourceSpace,optionalAlgoInput(hMatrixRegMat))(WeightingInSensorSpace,optionalAlgoInput(hMatrixNoiseCov))(matrixU,optionalAlgoInput(hMatrixU))(singularValues,optionalAlgoInput(hSingularValues))(matrixV,optionalAlgoInput(hMatrixV))) );
121120

‎src/Modules/Legacy/Inverse/SolveInverseProblemWithTikhonovTSVD.cc‎

Lines changed: 31 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -67,13 +67,21 @@ SolveInverseProblemWithTikhonovTSVD::SolveInverseProblemWithTikhonovTSVD() : Mod
6767

6868
void SolveInverseProblemWithTikhonovTSVD::setStateDefaults()
6969
{
70-
setStateStringFromAlgo(Parameters::TikhonovImplementation);
70+
auto state = get_state();
71+
72+
state->setValue( Parameters::TikhonovImplementation, std::string("TikhonovTSVD") );
73+
state->setValue(Parameters::LambdaMin,1);
74+
state->setValue(Parameters::LambdaMax,1);
75+
state->setValue(Parameters::LambdaNum, 1);
76+
state->setValue( Parameters::LambdaResolution, 1);
77+
78+
// setStateStringFromAlgo(Parameters::TikhonovImplementation);
7179
setStateStringFromAlgoOption(Parameters::RegularizationMethod);
7280
setStateDoubleFromAlgo(Parameters::LambdaFromDirectEntry);
73-
setStateDoubleFromAlgo(Parameters::LambdaMin);
74-
setStateDoubleFromAlgo(Parameters::LambdaMax);
75-
setStateIntFromAlgo(Parameters::LambdaNum);
76-
setStateDoubleFromAlgo(Parameters::LambdaResolution);
81+
// setStateDoubleFromAlgo(Parameters::LambdaMin);
82+
// setStateDoubleFromAlgo(Parameters::LambdaMax);
83+
// setStateIntFromAlgo(Parameters::LambdaNum);
84+
// setStateDoubleFromAlgo(Parameters::LambdaResolution);
7785
setStateDoubleFromAlgo(Parameters::LambdaSliderValue);
7886
setStateIntFromAlgo(Parameters::LambdaCorner);
7987
setStateStringFromAlgo(Parameters::LCurveText);
@@ -96,26 +104,39 @@ void SolveInverseProblemWithTikhonovTSVD::execute()
96104
auto hSingularValues = getOptionalInput(singularValues);
97105
auto hMatrixV = getOptionalInput(matrixV);
98106

99-
std::cout << "gato" << std::endl;
100107

101108
if (needToExecute())
102109
{
110+
111+
// Obtain rank of forward matrix
112+
int rank;
113+
if ( hSingularValues ) {
114+
rank = (*hSingularValues)->nrows();
115+
}
116+
else{
117+
Eigen::FullPivLU<SCIRun::Core::Datatypes::DenseMatrix::EigenBase> lu_decomp(*forward_matrix_h);
118+
rank = lu_decomp.rank();
119+
}
120+
121+
std::cout << "Computed rank: " << rank << std::endl;
103122
// set parameters
104123
auto state = get_state();
105-
// set parameters
106-
state->setValue( Parameters::TikhonovImplementation, std::string("TikhonovTSVD") );
124+
125+
// state->setValue( Parameters::TikhonovImplementation, std::string("TikhonovTSVD") );
107126
setAlgoStringFromState(Parameters::TikhonovImplementation);
108127
setAlgoOptionFromState(Parameters::RegularizationMethod);
109128
setAlgoDoubleFromState(Parameters::LambdaFromDirectEntry);
129+
// state->setValue(Parameters::LambdaMin,1);
130+
// state->setValue(Parameters::LambdaMax,double(rank)); // casting to double to keep consistency across tikhonov types
131+
// state->setValue(Parameters::LambdaNum, rank); // casting to double to keep consistency across tikhonov types
132+
// state->setValue( Parameters::LambdaResolution, 1);
110133
setAlgoDoubleFromState(Parameters::LambdaMin);
111134
setAlgoDoubleFromState(Parameters::LambdaMax);
112135
setAlgoIntFromState(Parameters::LambdaNum);
113-
setAlgoDoubleFromState(Parameters::LambdaResolution);
114136
setAlgoDoubleFromState(Parameters::LambdaSliderValue);
115137
setAlgoIntFromState(Parameters::LambdaCorner);
116138
setAlgoStringFromState(Parameters::LCurveText);
117139

118-
// run
119140
// run
120141
auto output = algo().run( withInputData((ForwardMatrix, forward_matrix_h)(MeasuredPotentials,hMatrixMeasDat)(MeasuredPotentials,hMatrixMeasDat)(WeightingInSourceSpace,optionalAlgoInput(hMatrixRegMat))(WeightingInSensorSpace,optionalAlgoInput(hMatrixNoiseCov))(matrixU,optionalAlgoInput(hMatrixU))(singularValues,optionalAlgoInput(hSingularValues))(matrixV,optionalAlgoInput(hMatrixV))) );
121142

0 commit comments

Comments
 (0)