Skip to content

Commit 84f0348

Browse files
jcollfontjcollfont
committed
Tikhonov seems to be fixed
1 parent 9c104f8 commit 84f0348

File tree

5 files changed

+69
-104
lines changed

5 files changed

+69
-104
lines changed

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

Lines changed: 0 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -70,7 +70,6 @@ using namespace SCIRun::Core::Algorithms::Inverse;
7070
//
7171
// A^-1 = M3 * G^-1 * M4
7272
//...........................................................................................................
73-
7473
const int sizeB = M1.ncols();
7574
const int sizeSolution = M3.nrows();
7675
const int numTimeSamples = y.ncols();
@@ -107,8 +106,6 @@ using namespace SCIRun::Core::Algorithms::Inverse;
107106
const int M = forwardMatrix_.nrows();
108107
const int N = forwardMatrix_.ncols();
109108

110-
int x = 0;
111-
112109
// PREALOCATE VARIABLES and MATRICES
113110
DenseMatrix forward_transpose = forwardMatrix_.transpose();
114111

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

Lines changed: 60 additions & 68 deletions
Original file line numberDiff line numberDiff line change
@@ -36,7 +36,7 @@ Last modification : April 20 2017
3636
#include <Core/Algorithms/Legacy/Inverse/TikhonovAlgoAbstractBase.h>
3737
#include <Core/Algorithms/Legacy/Inverse/TikhonovImpl.h>
3838
#include <Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithStandardTikhonovImpl.h>
39-
// #include <Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithTikhonovSVD_impl.h>
39+
#include <Core/Algorithms/Legacy/Inverse/SolveInverseProblemWithTikhonovSVD_impl.h>
4040

4141
// Datatypes
4242
#include <Core/Datatypes/Matrix.h>
@@ -87,14 +87,13 @@ TikhonovAlgoAbstractBase::TikhonovAlgoAbstractBase()
8787
{
8888
using namespace Parameters;
8989

90-
addParameter(TikhonovImplementation, "NoMethodSelected" );
91-
addOption(RegularizationMethod, "lcurve", "single|slider|lcurve");
92-
addParameter(regularizationChoice, 0);
93-
addParameter(LambdaFromDirectEntry,1e-6);
94-
// addParameter(lambdaDoubleSpinBox,1e-6);
90+
addParameter(Parameters::TikhonovImplementation, std::string("NoMethodSelected") );
91+
addOption(Parameters::RegularizationMethod, "lcurve", "single|slider|lcurve");
92+
addParameter(Parameters::regularizationChoice, 0);
93+
addParameter(Parameters::LambdaFromDirectEntry,1e-6);
9594
addParameter(Parameters::LambdaMin,1e-6);
9695
addParameter(Parameters::LambdaMax,1);
97-
addParameter(Parameters::LambdaNum,1000);
96+
addParameter(Parameters::LambdaNum,200);
9897
addParameter(Parameters::LambdaResolution,1e-6);
9998
addParameter(Parameters::LambdaSliderValue,0);
10099
addParameter(Parameters::LambdaCorner,0);
@@ -209,19 +208,25 @@ AlgorithmOutput TikhonovAlgoAbstractBase::run(const AlgorithmInput & input) cons
209208

210209
// Determine specific Tikhonov Implementation
211210
TikhonovImpl *algoImpl;
212-
std::cout << "Selecting Tikhonov: " << TikhonovImplementation_gotten << std::endl;
213-
if ( TikhonovImplementation_gotten == "standardTikhonov" ){
211+
if ( TikhonovImplementation_gotten == std::string("standardTikhonov") ){
212+
214213
// get Parameters
215214
int regularizationChoice_ = get(Parameters::regularizationChoice).toInt();
216215
int regularizationSolutionSubcase_ = get(Parameters::regularizationSolutionSubcase).toInt();
217216
int regularizationResidualSubcase_ = get(Parameters::regularizationResidualSubcase).toInt();
218217

219218
algoImpl = new SolveInverseProblemWithStandardTikhonovImpl( *castMatrix::toDense(forwardMatrix_), *castMatrix::toDense(measuredData_), *castMatrix::toDense(sourceWeighting_), *castMatrix::toDense(sensorWeighting_), regularizationChoice_, regularizationSolutionSubcase_, regularizationResidualSubcase_);
220219
}
221-
else if ( TikhonovImplementation_gotten == "TikhonovSVD" ){
222-
// algoImpl = new SolveInverseProblemWithTikhonovSVD_impl( *castMatrix::toDense(forwardMatrix_), *castMatrix::toDense(measuredData_), *castMatrix::toDense(sourceWeighting_), *castMatrix::toDense(sensorWeighting_));
220+
else if ( TikhonovImplementation_gotten == std::string("TikhonovSVD") ){
221+
222+
// get Parameters
223+
int regularizationChoice_ = get(Parameters::regularizationChoice).toInt();
224+
int regularizationSolutionSubcase_ = get(Parameters::regularizationSolutionSubcase).toInt();
225+
int regularizationResidualSubcase_ = get(Parameters::regularizationResidualSubcase).toInt();
226+
227+
algoImpl = new SolveInverseProblemWithTikhonovSVD_impl(*castMatrix::toDense(forwardMatrix_), *castMatrix::toDense(measuredData_), *castMatrix::toDense(sourceWeighting_), *castMatrix::toDense(sensorWeighting_));
223228
}
224-
else if ( TikhonovImplementation_gotten== "TikhonovTSVD" ){
229+
else if ( TikhonovImplementation_gotten == std::string("TikhonovTSVD") ){
225230
THROW_ALGORITHM_PROCESSING_ERROR("Tikhonov TSVD not implemented yet");
226231
}
227232
else{
@@ -247,6 +252,7 @@ AlgorithmOutput TikhonovAlgoAbstractBase::run(const AlgorithmInput & input) cons
247252
else if (RegularizationMethod_gotten == "lcurve")
248253
{
249254
lambda_ = computeLcurve( algoImpl, input );
255+
std::cout << "Lambda: " << lambda_ << std::endl;
250256
}
251257
else
252258
{
@@ -287,21 +293,18 @@ AlgorithmOutput TikhonovAlgoAbstractBase::run(const AlgorithmInput & input) cons
287293
/////// compute L-curve
288294
double TikhonovAlgoAbstractBase::computeLcurve( const SCIRun::Core::Algorithms::Inverse::TikhonovImpl * algoImpl, const AlgorithmInput & input ) const
289295
{
290-
291296
// get inputs
292297
auto forwardMatrix_ = input.get<Matrix>(TikhonovAlgoAbstractBase::ForwardMatrix);
293298
auto measuredData_ = input.get<Matrix>(TikhonovAlgoAbstractBase::MeasuredPotentials);
294299
auto sourceWeighting_ = input.get<Matrix>(TikhonovAlgoAbstractBase::WeightingInSourceSpace);
295300
auto sensorWeighting_ = input.get<Matrix>(TikhonovAlgoAbstractBase::WeightingInSensorSpace);
296-
297301
// define the step size of the lambda vector to be computed (distance between min and max divided by number of desired lambdas in log scale)
298302
const int nLambda = get(Parameters::LambdaNum).toInt();
299-
const int lambdaMin_ = get(Parameters::LambdaMin).toDouble();
300-
const int lambdaMax_ = get(Parameters::LambdaMax).toDouble();
301-
const double lam_step = pow(10.0, lambdaMax_ / lambdaMin_) / (nLambda-1);
303+
const double lambdaMin_ = get(Parameters::LambdaMin).toDouble();
304+
const double lambdaMax_ = get(Parameters::LambdaMax).toDouble();
305+
const double lam_step = (log10(lambdaMax_) - log10(lambdaMin_)) / (nLambda-1);
306+
std::cout << "Lambda power step: " << lam_step << ". Number: "<< nLambda <<". Lambda min: " << lambdaMin_ << ". Lambda max: "<< lambdaMax_<< ". Ratio: "<< lambdaMax_ / lambdaMin_ << std::endl;
302307
double lambda = 0;
303-
304-
const int sizeSolution = forwardMatrix_->ncols();
305308
double lambda_sq;
306309

307310
// prealocate vector of lambdas and eta and rho
@@ -314,15 +317,12 @@ double TikhonovAlgoAbstractBase::computeLcurve( const SCIRun::Core::Algorithms::
314317

315318
lambdaArray[0] = lambdaMin_;
316319

317-
// initialize counter
318-
int lambda_index = 0;
319-
320320
// for all lambdas
321321
for (int j = 0; j < nLambda; j++)
322322
{
323323
if (j)
324324
{
325-
lambdaArray[j] = lambdaArray[j-1] * lam_step;
325+
lambdaArray[j] = lambdaArray[j-1] * pow(10.0,lam_step);
326326
}
327327

328328
// set current lambda
@@ -377,13 +377,9 @@ double TikhonovAlgoAbstractBase::computeLcurve( const SCIRun::Core::Algorithms::
377377

378378
}
379379

380-
// // update L-curve
381-
// boost::shared_ptr<TikhonovAlgorithm::LCurveInput> lcurveInput(new TikhonovAlgorithm::LCurveInput(rho, eta, lambdaArray, nLambda));
382-
// lcurveInput_handle_ = lcurveInput;
383-
//
384-
// // Find corner in L-curve
385-
// lambda = FindCorner(*lcurveInput_handle_, lambda_index);
386-
//
380+
// Find corner in L-curve
381+
lambda = FindCorner( rho, eta, lambdaArray, nLambda );
382+
387383
// // update GUI
388384
// if (updateLCurveGui_)
389385
// updateLCurveGui_(lambda, *lcurveInput_handle_, lambda_index);
@@ -397,13 +393,9 @@ double TikhonovAlgoAbstractBase::computeLcurve( const SCIRun::Core::Algorithms::
397393

398394

399395
///// Find Corner, find the maximal curvature which corresponds to the L-curve corner
400-
double TikhonovAlgoAbstractBase::FindCorner( LCurveInput & Linput, const AlgorithmInput & input, int& lambda_index)
396+
double TikhonovAlgoAbstractBase::FindCorner( const std::vector<double>& rho, const std::vector<double>& eta, const std::vector<double>& lambdaArray, const int nLambda )
401397
{
402-
const int nLambda = Linput.nLambda_;
403-
const std::vector<double>& rho = Linput.rho_;
404-
const std::vector<double>& eta = Linput.eta_;
405-
const std::vector<double>& lambdaArray = Linput.lambdaArray_;
406-
398+
int lambda_index;
407399
std::vector<double> deta(nLambda);
408400
std::vector<double> ddeta(nLambda);
409401
std::vector<double> drho(nLambda);
@@ -450,38 +442,38 @@ double TikhonovAlgoAbstractBase::FindCorner( LCurveInput & Linput, const Algorit
450442
return lambdaArray[lambda_index];
451443
}
452444

453-
///// Search for closest Lambda to given lambda
454-
double TikhonovAlgoAbstractBase::LambdaLookup( LCurveInput& Linput, double lambda, int& lambda_index, const double epsilon)
455-
{
456-
const int nLambda = Linput.nLambda_;
457-
const std::vector<double>& lambdaArray = Linput.lambdaArray_;
458-
459-
for (int i = 0; i < nLambda-1; ++i)
460-
{
461-
if (i > 0 && (lambda < lambdaArray[i-1] || lambda > lambdaArray[i+1])) continue;
462-
463-
double lambda_step_midpoint = std::abs(lambdaArray[i+1] - lambdaArray[i])/2;
464-
465-
if (std::abs(lambda - lambdaArray[i]) <= epsilon) // TODO: is this a reasonable comparison???
466-
{
467-
lambda_index = i;
468-
return lambdaArray[lambda_index];
469-
}
470-
471-
if (std::abs(lambda - lambdaArray[i]) < lambda_step_midpoint)
472-
{
473-
lambda_index = i;
474-
return lambdaArray[lambda_index];
475-
}
476-
477-
if (std::abs(lambda - lambdaArray[i+1]) < lambda_step_midpoint)
478-
{
479-
lambda_index = i+1;
480-
return lambdaArray[lambda_index];
481-
}
482-
}
483-
return -1;
484-
}
445+
// ///// Search for closest Lambda to given lambda
446+
// double TikhonovAlgoAbstractBase::LambdaLookup( LCurveInput& Linput, double lambda, int& lambda_index, const double epsilon)
447+
// {
448+
// const int nLambda = Linput.nLambda_;
449+
// const std::vector<double>& lambdaArray = Linput.lambdaArray_;
450+
//
451+
// for (int i = 0; i < nLambda-1; ++i)
452+
// {
453+
// if (i > 0 && (lambda < lambdaArray[i-1] || lambda > lambdaArray[i+1])) continue;
454+
//
455+
// double lambda_step_midpoint = std::abs(lambdaArray[i+1] - lambdaArray[i])/2;
456+
//
457+
// if (std::abs(lambda - lambdaArray[i]) <= epsilon) // TODO: is this a reasonable comparison???
458+
// {
459+
// lambda_index = i;
460+
// return lambdaArray[lambda_index];
461+
// }
462+
//
463+
// if (std::abs(lambda - lambdaArray[i]) < lambda_step_midpoint)
464+
// {
465+
// lambda_index = i;
466+
// return lambdaArray[lambda_index];
467+
// }
468+
//
469+
// if (std::abs(lambda - lambdaArray[i+1]) < lambda_step_midpoint)
470+
// {
471+
// lambda_index = i+1;
472+
// return lambdaArray[lambda_index];
473+
// }
474+
// }
475+
// return -1;
476+
// }
485477

486478
// ////////// update L-curve graph
487479
// void TikhonovAlgoAbstractBase::update_graph( LCurveInput & input, double lambda, int lambda_index, const double epsilon)

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

Lines changed: 2 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -71,12 +71,6 @@ namespace Inverse {
7171
static const AlgorithmOutputName RegInverse;
7272

7373
// Define algorithm choices
74-
enum TikhonovImplementationChoice {
75-
NoMethodSelected,
76-
standardTikhonov,
77-
TikhonovSVD,
78-
TikhonovTSVD
79-
};
8074
enum AlgorithmChoice {
8175
automatic,
8276
underdetermined,
@@ -91,17 +85,6 @@ namespace Inverse {
9185
residual_constrained_squared
9286
};
9387

94-
struct SCISHARE LCurveInput
95-
{
96-
const std::vector<double> rho_;
97-
const std::vector<double> eta_;
98-
const std::vector<double> lambdaArray_;
99-
const int nLambda_;
100-
101-
LCurveInput(const std::vector<double>& rho, const std::vector<double>& eta, const std::vector<double>& lambdaArray, const int nLambda);
102-
};
103-
104-
10588
// constructor
10689
TikhonovAlgoAbstractBase();
10790

@@ -110,8 +93,8 @@ namespace Inverse {
11093

11194
// defined public functions
11295
// void update_graph( const AlgorithmInput & input, double lambda, int lambda_index, double epsilon);
113-
static double FindCorner( LCurveInput & Linput, const AlgorithmInput & input, int& lambda_index);
114-
static double LambdaLookup(LCurveInput& input, double lambda, int& lambda_index, const double epsilon);
96+
static double FindCorner( const std::vector<double>& rho, const std::vector<double>& eta, const std::vector<double>& lambdaArray, const int nLambda );
97+
// static double LambdaLookup(LCurveInput& input, double lambda, int& lambda_index, const double epsilon);
11598
double computeLcurve( const SCIRun::Core::Algorithms::Inverse::TikhonovImpl * algoImpl, const AlgorithmInput & input ) const;
11699

117100
bool checkInputMatrixSizes( const AlgorithmInput & input ) const;

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

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -43,7 +43,7 @@ SolveInverseProblemWithTikhonovDialog::SolveInverseProblemWithTikhonovDialog(con
4343
setupUi(this);
4444
setWindowTitle(QString::fromStdString(name));
4545
fixSize();
46-
46+
4747
GuiStringTranslationMap lambdaMethod_;
4848
lambdaMethod_.insert(StringPair("Direct entry", "single"));
4949
lambdaMethod_.insert(StringPair("Slider", "slider"));
@@ -68,7 +68,7 @@ SolveInverseProblemWithTikhonovDialog::SolveInverseProblemWithTikhonovDialog(con
6868

6969
addComboBoxManager(lambdaMethodComboBox_, Parameters::RegularizationMethod, lambdaMethod_);
7070
addTextEditManager(lCurveTextEdit_, Parameters::LCurveText);
71-
71+
7272
connect(lambdaSlider_, SIGNAL(valueChanged(int)), this, SLOT(setSpinBoxValue(int)));
7373
connect(lambdaSliderDoubleSpinBox_, SIGNAL(valueChanged(double)), this, SLOT(setSliderValue(double)));
7474
connect(lambdaMinDoubleSpinBox_, SIGNAL(valueChanged(double)), this, SLOT(setSliderMin(double)));

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

Lines changed: 5 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -67,7 +67,7 @@ void SolveInverseProblemWithTikhonov::setStateDefaults()
6767
setStateDoubleFromAlgo(Parameters::LambdaFromDirectEntry);
6868
setStateDoubleFromAlgo(Parameters::LambdaMin);
6969
setStateDoubleFromAlgo(Parameters::LambdaMax);
70-
setStateDoubleFromAlgo(Parameters::LambdaNum);
70+
setStateIntFromAlgo(Parameters::LambdaNum);
7171
setStateDoubleFromAlgo(Parameters::LambdaResolution);
7272
setStateDoubleFromAlgo(Parameters::LambdaSliderValue);
7373
setStateIntFromAlgo(Parameters::LambdaCorner);
@@ -91,32 +91,25 @@ void SolveInverseProblemWithTikhonov::execute()
9191
if (needToExecute())
9292
{
9393

94-
9594
auto state = get_state();
9695
// set parameters
97-
std::cout << "The parameter was set to: " << state->getValue(Parameters::TikhonovImplementation).toString() << std::endl;
98-
std::cout << "Selecting Tikhonov in the module: " ;
99-
100-
state->setValue( Parameters::TikhonovImplementation, "standardTikhonov" );
96+
state->setValue( Parameters::TikhonovImplementation, std::string("standardTikhonov") );
97+
setAlgoStringFromState(Parameters::TikhonovImplementation);
10198
setAlgoOptionFromState(Parameters::RegularizationMethod);
10299
setAlgoIntFromState(Parameters::regularizationChoice);
103100
setAlgoDoubleFromState(Parameters::LambdaFromDirectEntry);
104101
setAlgoDoubleFromState(Parameters::LambdaMin);
105102
setAlgoDoubleFromState(Parameters::LambdaMax);
106-
setAlgoDoubleFromState(Parameters::LambdaNum);
103+
setAlgoIntFromState(Parameters::LambdaNum);
107104
setAlgoDoubleFromState(Parameters::LambdaResolution);
108105
setAlgoDoubleFromState(Parameters::LambdaSliderValue);
109106
setAlgoIntFromState(Parameters::LambdaCorner);
110107
setAlgoStringFromState(Parameters::LCurveText);
111108
setAlgoIntFromState(Parameters::regularizationSolutionSubcase);
112109
setAlgoIntFromState(Parameters::regularizationResidualSubcase);
113110

114-
115-
std::cout << state->getValue(Parameters::TikhonovImplementation).toString() << std::endl;
116-
117111
// run
118-
// auto output = algo().run( withInputData((ForwardMatrix, forward_matrix_h)(MeasuredPotentials,hMatrixMeasDat)(WeightingInSourceSpace,hMatrixRegMat)(WeightingInSensorSpace, hMatrixNoiseCov)));
119-
auto output = algo().run( withInputData((ForwardMatrix, forward_matrix_h)(MeasuredPotentials,hMatrixMeasDat)) );
112+
auto output = algo().run( withInputData((ForwardMatrix, forward_matrix_h)(MeasuredPotentials,hMatrixMeasDat)(MeasuredPotentials,hMatrixMeasDat)(WeightingInSourceSpace,optionalAlgoInput(hMatrixRegMat))(WeightingInSensorSpace,optionalAlgoInput(hMatrixNoiseCov))) );
120113

121114
// update L-curve
122115
/* NO EXISTE

0 commit comments

Comments
 (0)