[Insight-developers] (no subject)
C. Aaron Cois
cacst11+ at pitt . edu
Thu, 20 Jun 2002 14:37:00 -0400
This is a multi-part message in MIME format.
------=_NextPart_000_001A_01C21867.EFCB9570
Content-Type: text/plain;
charset="iso-8859-1"
Content-Transfer-Encoding: 7bit
Lydia, you are right. I apologize, and here is our
itkLevenbergMarquardtOptimizerTest.cxx.
----- Original Message -----
From: "Lydia Ng" <lng@insightful.com>
To: "C. Aaron Cois" <cacst11@pitt.edu>;
<insight-developers@public.kitware.com>
Sent: Thursday, June 20, 2002 1:32 PM
Subject: RE: [Insight-developers] (no subject)
Hi Aaron,
The error message indicates that something has written
to memory off the end of some allocation.
Could you post your
itkLevenbergMarquardtOptimizerTest.cxx
I think you sent out
itkLevenbergMarquardtOptimizer.cxx
in error.
- Lydia
-----Original Message-----
From: C. Aaron Cois [mailto:cacst11+@pitt.edu]
Sent: Thursday, June 20, 2002 10:28 AM
To: insight-developers@public.kitware.com
Subject: [Insight-developers] (no subject)
We have run across an error that has stumped us for almost a week. It seems
like a good time to see if any other developers have any thoughts or similar
experiences.
We modified itkLevenbergMarquardtOptimizerTest.cxx to operate on our
particular cost function. The major modifications were:
1. GetValue() modified
2. GetDerivative() modified
3. Added ERF function
4. Changed Max_Iterations to 100
And a few other changes like prints and const. We used data generated with
known parameters from excel and the test will fit the parameters almost
perfectly:
----------------
vnl_levenberg_marquardt: 13 iterations, 66 evaluations, 5 residuals. RMS
error
start/end 60.749/0.000208907
End condition = Failed X Tolerance too small
Number of iters = 13
Number of evals = 66
Solution = (31.9992,255.022,4.0002,2.00008)
Test passed.
----------------
After this solutions is printed out, we get a crash as the smart pointer for
the cost function goes out of scope...some kind of crash on cleanup.
I've included the source file (replace
itkLevenbergMarquardtOptimizerTest.cxx and hit 4 when running) and a screen
capture of the error.
Any ideas?
Thanks
_______________________________________________
Insight-developers mailing list
Insight-developers@public.kitware.com
http://public.kitware.com/mailman/listinfo/insight-developers
------=_NextPart_000_001A_01C21867.EFCB9570
Content-Type: application/octet-stream;
name="itkLevenbergMarquardtOptimizerTest.cxx"
Content-Transfer-Encoding: quoted-printable
Content-Disposition: attachment;
filename="itkLevenbergMarquardtOptimizerTest.cxx"
/*=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: itkLevenbergMarquardtOptimizerTest.cxx,v $
Language: C++
Date: $Date: 2002/03/07 21:04:15 $
Version: $Revision: 1.20 $
Copyright (c) 2002 Insight Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for =
details.
This software is distributed WITHOUT ANY WARRANTY; without even=20
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR =
PURPOSE. See the above copyright notices for more information.
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=
=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D=3D*/
#include <itkLevenbergMarquardtOptimizer.h>
#include <vnl/vnl_vector.h>
#include <vnl/vnl_matrix.h>
#include <vnl/vnl_math.h>
#include <itkArray.h>
#include <itkArray2D.h>
const double INV_SQRT_TWO_PI =3D .398942280401;
const double SQUARE_ROOT_OF_TWO =3D 1.41421356237;
typedef vnl_matrix<double> MatrixType;
typedef vnl_vector<double> VectorType;
/**=20
*
* This example minimize the equation:
*
* ( a * x + b * y + c )=20
* -( 2 * x + 3 * y + 4 )
* =20
* for the (a,b,c) parameters
*
* the solution is the vector | 2 3 4 |
*
* (x,y) values are sampled over a rectangular =
doitkLevenbergMarquardtOptimizerTest
* whose size is defined by XRange and YRange
*
*/=20
class LMCostFunction : public itk::MultipleValuedCostFunction
{
public:
typedef LMCostFunction Self;
typedef itk::MultipleValuedCostFunction Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
itkNewMacro( Self );
enum { XRange =3D 1,
YRange =3D 1 }; // size of the =
doitkLevenbergMarquardtOptimizerTest to sample the cost function
=20
enum { SpaceDimension =3D 4 };
enum { RangeDimension =3D 5 };
typedef Superclass::ParametersType ParametersType;
typedef Superclass::DerivativeType DerivativeType;
typedef Superclass::MeasureType MeasureType;
LMCostFunction():
m_Measure(RangeDimension),
m_Derivative(SpaceDimension,RangeDimension),
m_TheoreticalData(SpaceDimension) =20
{
m_Measure.resize(RangeDimension);
m_Derivative.resize(SpaceDimension,RangeDimension);
m_TheoreticalData.resize(RangeDimension);
=20
m_TheoreticalData[0] =3D 37.80128449;
m_TheoreticalData[1] =3D 49.03584331;
m_TheoreticalData[2] =3D 72.45708986;
m_TheoreticalData[3] =3D 110.6770725;
m_TheoreticalData[4] =3D 159.5;
=20
std::cout << std::endl;
}
MeasureType GetValue( const ParametersType & parameters ) const
{ =20
double a =3D parameters[0]; // lower intensity
double b =3D parameters[1]; // upper intensity
double c =3D parameters[2]; // mean
double d =3D parameters[3]; // standard deviation
// Compute points of the function over a square domain
unsigned valueindex =3D 0;
double arg =3D 0;
for( int x =3D 0; x<=3DRangeDimension; x++ )=20
{
const double xd =3D (double)x;
=20
arg =3D ((xd-c)/(d*SQUARE_ROOT_OF_TWO));
m_Measure[valueindex] =3D (a+.5*(b)*(1 + erf(arg)));
m_Measure[valueindex] -=3D m_TheoreticalData[valueindex];
std::cout << m_Measure[valueindex] << " ";
valueindex++;
}
std::cout << std::endl;
return m_Measure;=20
}
void GetDerivative( const ParametersType & parameters,
DerivativeType & derivative ) const
{
=20
double a =3D parameters[0];
double b =3D parameters[1];
double c =3D parameters[2];
double d =3D parameters[3];
// Compute points of the function over a square domain
unsigned valueindex =3D 0;
=20
double arg =3D 0;
for( int x =3D 0; x<=3DRangeDimension; x++ )=20
{
const double xd =3D (double)x;
arg =3D ((x-c)/(d*SQUARE_ROOT_OF_TWO));
m_Derivative[0][valueindex] =3D 1;
m_Derivative[1][valueindex] =3D (.5*(1 + erf(arg)));
m_Derivative[2][valueindex] =3D =
-(INV_SQRT_TWO_PI*b*(1/d)*(exp(arg*arg)));
m_Derivative[3][valueindex] =3D -(INV_SQRT_TWO_PI*b*(x - =
c)*(1/(d*d))*exp(arg*arg));
valueindex++;
}
for(unsigned int dim1=3D0; dim1 < SpaceDimension; dim1++)
{
std::cout << m_Derivative[dim1] << " ";
}
std::cout << std::endl;
}
unsigned int GetNumberOfParameters(void) const
{
return SpaceDimension;
}
unsigned int GetNumberOfValues(void) const
{
return RangeDimension;
}
double erf(double argument) const
{
int temp =3D 0;
double erfValue =3D 0;
double slope;
=20
double y[300] =3D=20
{
0, .011283416, .022564575, .033841222, .045111106, =
.056371978, .067621594, .07885772, .090078126, .101280594,
.112462916, .123622896, .134758352, .145867115, .156947033, =
.167995971, .179011813, .189992461, .200935839, .211839892,
.222702589, .233521923, .244295911, .255022599, .265700058, =
.276326389, .286899723, .297418219, .307880068, .318283496,
.328626759, .33890815, .349125995, .359278655, .369364529, =
.379382053, .3893297, .399205983, .409009452, .418738698,
.428392352, .43796909, .447467618, .456886694, .466225115, =
.475481719, .484655389, .49374505, .50274967, .51166826,
.520499876, .529243617, .537898627, .546464093, .554939245, =
.563323359, .571615763, .579815806, .5879229, .595936496,
.60385609, .611681217, .61941146, .627046441, .634585826, =
.642029324, .649376683, .656627696, .663782195, .670840052,
.677801193, .684665264, .691432825, .698103704, .704677825, =
.71115543, .717536534, .723821437, .730010238, .73610324,
.74210079, .748003138, .75381059, .759523625, .76514256, =
.770667933, .776100122, .781439725, .786687219, .791843127,
.796908113, .801882743, .80676762, .811563474, .816270948, =
.820890718, .825423575, .82987023, .834231422, .838508001,
.842700735, .846810448, .850837952, .854784156, .8586499, =
.862436067, .866143531, .86977325, .873326119, .876803068,
.880205041, .88353297, .886787854, .88997064, .893082302, =
.896123821, .899096169, .90200037, .904837402, .907608265,
.91031396, .912955492, .915533856, .918050082, .920505165, =
.922900112, .925235928, .927513617, .929734183, .931898615,
.934007929, .936063109, .938065143, .940015016, .941913707, =
.943762189, .94556143, .947312386, .949016025, .950673287,
.952285112, .953852432, .955376173, .956857248, .958296565, =
.959695022, .961053506, .962372893, .963654059, .964897859,
.966105142, .967276744, .968413493, .969516206, .970585687, =
.971622731, .97262812, .973602626, .974547008, .975462012,
.97634838, .977206834, .978038086, .978842837, .979621778, =
.980375583, .98110492, .98181044, .982492786, .983152586,
.983790458, .984407007, .985002827, .985578499, .986134593, =
.98667167, .987190274, .987690941, .988174195, .988640548,
.989090501, .989524544, .989943156, .990346805, .990735947, =
.99111103, .991472488, .991820747, .992156222, .992479318,
.992790429, .99308994, .993378225, .99365565, .99392257, =
.994179333, .994426275, .994663724, .994892, .995111413,
.995322265, .995524849, .995719451, .995906348, .996085809, =
.996258096, .996423462, .996582153, .996734409, .99688046,
.997020533, .997154845, .997283607, .997407023, .997525293, =
.997638607, .997747152, .997851108, .997950649, .998045943,
.998137154, .998224438, .998307948, .998387832, .998464231, =
.998537283, .998607121, .998673872, .998737661, .998798606, =
=20
.998856823, .998912423, .998965513, .999016195, .99906457, =
.999110733, .999154777, .99919679, .999236858, .999275064,
.999311486, .999346202, .999379283, .999410802, .999440826, =
.99946942, .999496646, .999522566, .999547236, .999570712,
.999593048, .999614295, .999634501, .999653714, .999671979, =
.99968934, .999705837, .999721511, .9997364, .999750539,
.999763966, .999776711, .999788809, .999800289, .999811181, =
.999821512, .999831311, .999840601, .999849409, .999857757,
.999865667, .999873162, .999880261, .999886985, .999893351, =
.999899378, .999905082, .99991048, .999915587, .999920418,
.999924987, .999929307, .99993339, .99993725, .999940898, =
.999944344, .999947599, .999950673, .999953576, .999956316,
.999958902, .999961343, .999963645, .999965817, .999967866, =
.999969797, .999971618, .999973334, .999974951, .999976474
};
=20
if(argument < -3 || argument > 3)
{
if(argument > 0)
{
erfValue =3D 1;
}
else
{
erfValue =3D -1;
}
}
// interpolation between table lookup entries
else
{
if(argument > 0)
{ =20
temp =3D (int)(argument * 100); =20
if(argument =3D=3D (int)temp) erfValue =3D .999976474;
else
{
slope =3D (y[temp + 1] - y[temp])/(((float)temp + 1)/100 - =
((float)temp/100));
erfValue =3D slope * (argument - ((float)temp + 1)/100) + =
y[temp+1]; =20
}
}=20
else
{
temp =3D -(int)(argument * 100);
slope =3D (-y[temp + 1] + y[temp])/(-((float)temp + 1)/100 + =
((float)temp/100));
erfValue =3D (slope * (argument + ((float)temp + 1)/100) - =
y[temp+1]);
}
}
return erfValue;
}
private:
mutable MeasureType m_Measure;
mutable DerivativeType m_Derivative;
MeasureType m_TheoreticalData;
};
int itkLevenbergMarquardtOptimizerTest(int, char**)=20
{
std::cout << "Levenberg Marquardt optimizer test \n \n";=20
=20
typedef itk::LevenbergMarquardtOptimizer OptimizerType;
typedef OptimizerType::InternalOptimizerType vnlOptimizerType;
=20
=20
// Declaration of a itkOptimizer
OptimizerType::Pointer Optimizer =3D OptimizerType::New();
// Declaration of the CostFunction adaptor
LMCostFunction::Pointer costFunction =3D LMCostFunction::New();
typedef LMCostFunction::ParametersType ParametersType;
ParametersType parameters(LMCostFunction::SpaceDimension);
costFunction->GetValue(parameters);
=20
std::cout << costFunction->GetNumberOfValues() << "\n";
try=20
{
Optimizer->SetCostFunction( costFunction.GetPointer() );
}
catch( itk::ExceptionObject & e )
{
std::cout << "Exception thrown ! " << std::endl;
std::cout << "An error ocurred during Optimization" << std::endl;
std::cout << e << std::endl;
return EXIT_FAILURE;
}
=20
const double F_Tolerance =3D 1e-15; // Function value tolerance
const double G_Tolerance =3D 1e-17; // Gradient magnitude =
tolerance=20
const double X_Tolerance =3D 1e-16; // Search space tolerance
const double Epsilon_Function =3D 1e-10; // Step
const int Max_Iterations =3D 100; // Maximum number of =
iterations
vnlOptimizerType * vnlOptimizer =3D Optimizer->GetOptimizer();
vnlOptimizer->set_f_tolerance( F_Tolerance );
vnlOptimizer->set_g_tolerance( G_Tolerance );
vnlOptimizer->set_x_tolerance( X_Tolerance );=20
vnlOptimizer->set_epsilon_function( Epsilon_Function );
vnlOptimizer->set_max_function_evals( Max_Iterations );
// We start not so far from the solution=20
typedef LMCostFunction::ParametersType ParametersType;
ParametersType initialValue(LMCostFunction::SpaceDimension);
initialValue[0] =3D 16;
initialValue[1] =3D 196;
initialValue[2] =3D 1;
initialValue[3] =3D 2;
OptimizerType::ParametersType =
currentValue(LMCostFunction::SpaceDimension);
currentValue =3D initialValue;
Optimizer->SetInitialPosition( currentValue );
try=20
{
Optimizer->StartOptimization();
}
catch( itk::ExceptionObject & e )
{
std::cout << "Exception thrown ! " << std::endl;
std::cout << "An error ocurred during Optimization" << std::endl;
std::cout << "Location =3D " << e.GetLocation() << std::endl;
std::cout << "Description =3D " << e.GetDescription() << std::endl;
return EXIT_FAILURE;
}
// Error codes taken from vxl/vnl/vnl_nonlinear_minimizer.h
std::cout << "End condition =3D ";
switch( vnlOptimizer->get_failure_code() )
{
case vnl_nonlinear_minimizer::ERROR_FAILURE:=20
std::cout << " Error Failure"; break;
case vnl_nonlinear_minimizer::ERROR_DODGY_INPUT:=20
std::cout << " Error Dogy Input"; break;
case vnl_nonlinear_minimizer::CONVERGED_FTOL:=20
std::cout << " Converged F Tolerance"; break;
case vnl_nonlinear_minimizer::CONVERGED_XTOL:=20
std::cout << " Converged X Tolerance"; break;
case vnl_nonlinear_minimizer::CONVERGED_XFTOL:
std::cout << " Converged XF Tolerance"; break;
case vnl_nonlinear_minimizer::CONVERGED_GTOL:=20
std::cout << " Converged G Tolerance"; break;
case vnl_nonlinear_minimizer::FAILED_TOO_MANY_ITERATIONS:
std::cout << " Too many iterations "; break;
case vnl_nonlinear_minimizer::FAILED_FTOL_TOO_SMALL:
std::cout << " Failed F Tolerance too small "; =
break;
case vnl_nonlinear_minimizer::FAILED_XTOL_TOO_SMALL:
std::cout << " Failed X Tolerance too small "; =
break;
case vnl_nonlinear_minimizer::FAILED_GTOL_TOO_SMALL:
std::cout << " Failed G Tolerance too small "; =
break;
}
std::cout << std::endl;
std::cout << "Number of iters =3D " << =
vnlOptimizer->get_num_iterations() << std::endl;
std::cout << "Number of evals =3D " << =
vnlOptimizer->get_num_evaluations() << std::endl; =20
std::cout << std::endl;
OptimizerType::ParametersType finalPosition;
finalPosition =3D Optimizer->GetCurrentPosition();
std::cout << "Solution =3D (";
std::cout << finalPosition[0] << "," ;
std::cout << finalPosition[1] << "," ;
std::cout << finalPosition[2] << "," ;
std::cout << finalPosition[3] << ")" << std::endl; =20
std::cout << "Test passed." << std::endl;
return EXIT_SUCCESS;
}
------=_NextPart_000_001A_01C21867.EFCB9570--