Skip to content

Commit

Permalink
ENH: Eliminates variable construct/destruct - important for large par…
Browse files Browse the repository at this point in the history 
…ameters spaces, e.g., during bspline transform optimization
  • Loading branch information
@aylward
aylward committed Apr 23, 2008
1 parent 778c5c0 commit a7b0a73
Show file tree
Hide file tree
Showing 5 changed files with 103 additions and 41 deletions.
38 changes: 23 additions & 15 deletions Code/Numerics/itkFRPROptimizer.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,15 @@ ::GetValueAndDerivative(ParametersType & p, double * val,
void
FRPROptimizer
::LineOptimize(ParametersType * p, ParametersType & xi, double * val)
{
ParametersType tempCoord( this->GetSpaceDimension() );
this->LineOptimize( p, xi, val, tempCoord );
}

void
FRPROptimizer
::LineOptimize(ParametersType * p, ParametersType & xi, double * val,
ParametersType & tempCoord )
{
this->SetLine(*p, xi);

Expand All @@ -79,15 +88,14 @@ ::LineOptimize(ParametersType * p, ParametersType & xi, double * val)
double bx;
double fb;

ParametersType pp = (*p);

this->LineBracket(&ax, &xx, &bx, &fa, &fx, &fb);
this->LineBracket(&ax, &xx, &bx, &fa, &fx, &fb, tempCoord);
this->SetCurrentLinePoint(xx, fx);

double extX = 0;
double extVal = 0;

this->BracketedLineOptimize(ax, xx, bx, fa, fx, fb, &extX, &extVal);
this->BracketedLineOptimize(ax, xx, bx, fa, fx, fb, &extX, &extVal,
tempCoord);
this->SetCurrentLinePoint(extX, extVal);

(*p) = this->GetCurrentPosition();
Expand All @@ -110,21 +118,21 @@ ::StartOptimization()

this->SetSpaceDimension(m_CostFunction->GetNumberOfParameters());

const unsigned int SpaceDimension = this->GetSpaceDimension();
FRPROptimizer::ParametersType tempCoord( this->GetSpaceDimension() );

double gg, gam, dgg;
FRPROptimizer::ParametersType g( SpaceDimension );
FRPROptimizer::ParametersType h( SpaceDimension );
FRPROptimizer::ParametersType xi( SpaceDimension );
FRPROptimizer::ParametersType g( this->GetSpaceDimension() );
FRPROptimizer::ParametersType h( this->GetSpaceDimension() );
FRPROptimizer::ParametersType xi( this->GetSpaceDimension() );

FRPROptimizer::ParametersType p( SpaceDimension );
FRPROptimizer::ParametersType p( this->GetSpaceDimension() );
p = this->GetInitialPosition();
this->SetCurrentPosition(p);

double fp;
this->GetValueAndDerivative(p, &fp, &xi);

for( i = 0; i < SpaceDimension; i++ )
for( i = 0; i < this->GetSpaceDimension(); i++ )
{
g[i] = -xi[i];
xi[i] = g[i];
Expand All @@ -141,12 +149,12 @@ ::StartOptimization()

double fret;
fret = fp;
this->LineOptimize(&p, xi, &fret);
this->LineOptimize(&p, xi, &fret, tempCoord);

if ( 2.0 * vcl_abs(fret - fp) <=
this->GetValueTolerance() * (vcl_abs(fret)+ vcl_abs(fp) + FRPR_TINY) )
{
if( limitCount < SpaceDimension )
if( limitCount < this->GetSpaceDimension() )
{
this->GetValueAndDerivative(p, &fp, &xi);
xi[limitCount] = 1;
Expand All @@ -170,15 +178,15 @@ ::StartOptimization()

if( m_OptimizationType == PolakRibiere )
{
for( i=0; i< SpaceDimension; i++ )
for( i=0; i< this->GetSpaceDimension(); i++ )
{
gg += g[i] * g[i];
dgg += (xi[i] + g[i]) * xi[i];
}
}
if( m_OptimizationType == FletchReeves )
{
for( i=0; i< SpaceDimension; i++ )
for( i=0; i< this->GetSpaceDimension(); i++ )
{
gg += g[i] * g[i];
dgg += xi[i] * xi[i];
Expand All @@ -193,7 +201,7 @@ ::StartOptimization()
}

gam = dgg/gg;
for( i = 0; i < SpaceDimension; i++)
for( i = 0; i < this->GetSpaceDimension(); i++)
{
g[i] = -xi[i];
xi[i] = g[i] + gam * h[i];
Expand Down
3 changes: 3 additions & 0 deletions Code/Numerics/itkFRPROptimizer.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,6 +96,9 @@ class ITK_EXPORT FRPROptimizer:

virtual void LineOptimize(ParametersType * p, ParametersType & xi,
double * val );
virtual void LineOptimize(ParametersType * p, ParametersType & xi,
double * val,
ParametersType & tempCoord );


private:
Expand Down
6 changes: 5 additions & 1 deletion Code/Numerics/itkOptimizer.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,12 +89,16 @@ class ITK_EXPORT Optimizer : public Object

bool m_ScalesInitialized;

// Keep m_CurrentPosition as a protected var so that subclasses can
// have fast access. This is important when optimizing high-dimensional
// spaces, e.g. bspline transforms.
ParametersType m_CurrentPosition;

private:
Optimizer(const Self&); //purposely not implemented
void operator=(const Self&); //purposely not implemented

ParametersType m_InitialPosition;
ParametersType m_CurrentPosition;
ScalesType m_Scales;

};
Expand Down
77 changes: 53 additions & 24 deletions Code/Numerics/itkPowellOptimizer.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,43 +55,47 @@ PowellOptimizer
::SetLine(const PowellOptimizer::ParametersType & origin,
const vnl_vector<double> & direction)
{
m_LineOrigin = origin;
m_LineDirection = direction;
for(unsigned int i=0; i<m_SpaceDimension; i++)
{
m_LineDirection[i] = m_LineDirection[i] / this->GetScales()[i];
m_LineOrigin[i] = origin[i];
m_LineDirection[i] = direction[i] / this->GetScales()[i];
}
}

double
PowellOptimizer
::GetLineValue(double x) const
{
PowellOptimizer::ParametersType xCoord( m_SpaceDimension );
PowellOptimizer::ParametersType tempCoord( m_SpaceDimension );
return this->GetLineValue(x, tempCoord);
}

double
PowellOptimizer
::GetLineValue(double x, ParametersType & tempCoord) const
{
for(unsigned int i=0; i<m_SpaceDimension; i++)
{
xCoord[i] = this->m_LineOrigin[i] + x * this->m_LineDirection[i];
tempCoord[i] = this->m_LineOrigin[i] + x * this->m_LineDirection[i];
}
if(m_Maximize)
{
return -(this->m_CostFunction->GetValue(xCoord));
return -(this->m_CostFunction->GetValue(tempCoord));
}
else
{
return this->m_CostFunction->GetValue(xCoord);
return this->m_CostFunction->GetValue(tempCoord);
}
}

void
PowellOptimizer
::SetCurrentLinePoint(double x, double fx)
{
PowellOptimizer::ParametersType xCoord( m_SpaceDimension );
for(unsigned int i=0; i<m_SpaceDimension; i++)
{
xCoord[i] = this->m_LineOrigin[i] + x * this->m_LineDirection[i];
this->m_CurrentPosition[i] = this->m_LineOrigin[i] + x * this->m_LineDirection[i];
}
this->SetCurrentPosition(xCoord);
if(m_Maximize)
{
this->SetCurrentCost(-fx);
Expand All @@ -100,6 +104,7 @@ ::SetCurrentLinePoint(double x, double fx)
{
this->SetCurrentCost(fx);
}
this->Modified();
}

void
Expand Down Expand Up @@ -143,6 +148,16 @@ void
PowellOptimizer
::LineBracket(double * x1, double * x2, double * x3,
double * f1, double * f2, double * f3)
{
PowellOptimizer::ParametersType tempCoord( m_SpaceDimension );
this->LineBracket( x1, x2, x3, f1, f2, f3, tempCoord);
}

void
PowellOptimizer
::LineBracket(double * x1, double * x2, double * x3,
double * f1, double * f2, double * f3,
ParametersType & tempCoord)
{
//
// Compute the golden ratio as a constant to be
Expand All @@ -153,7 +168,7 @@ ::LineBracket(double * x1, double * x2, double * x3,
//
// Get the value of the function for point x2
//
*f2 = this->GetLineValue( *x2 );
*f2 = this->GetLineValue( *x2, tempCoord );

//
// Compute extrapolated point using the golden ratio
Expand All @@ -162,7 +177,7 @@ ::LineBracket(double * x1, double * x2, double * x3,
{
// compute x3 on the side of x2
*x3 = *x1 + goldenRatio * ( *x2 - *x1 );
*f3 = this->GetLineValue( *x3 );
*f3 = this->GetLineValue( *x3, tempCoord );

// If the new point is a minimum
// then continue extrapolating in
Expand All @@ -174,14 +189,14 @@ ::LineBracket(double * x1, double * x2, double * x3,
*x2 = *x3;
*f2 = *f3;
*x3 = *x1 + goldenRatio * ( *x2 - *x1 );
*f3 = this->GetLineValue( *x3 );
*f3 = this->GetLineValue( *x3, tempCoord );
}
}
else
{
// compute x3 on the side of x1
*x3 = *x2 + goldenRatio * ( *x1 - *x2 );
*f3 = this->GetLineValue( *x3 );
*f3 = this->GetLineValue( *x3, tempCoord );

// If the new point is a minimum
// then continue extrapolating in
Expand All @@ -193,7 +208,7 @@ ::LineBracket(double * x1, double * x2, double * x3,
*x1 = *x3;
*f1 = *f3;
*x3 = *x2 + goldenRatio * ( *x1 - *x2 );
*f3 = this->GetLineValue( *x3 );
*f3 = this->GetLineValue( *x3, tempCoord );
}
}

Expand Down Expand Up @@ -226,7 +241,7 @@ ::LineBracket(double * x1, double * x2, double * x3,
// is no need for using it as a factor.
//
*x4 = *x1 - *x2 + *x3;
*f4 = this->GetLineValue( *x4 );
*f4 = this->GetLineValue( *x4, tempCoord );

//
// If the new value f4 at x4 is larger than f2
Expand Down Expand Up @@ -280,8 +295,19 @@ ::LineBracket(double * x1, double * x2, double * x3,
void
PowellOptimizer
::BracketedLineOptimize(double ax, double bx, double cx,
double itkNotUsed(fa), double functionValueOfb, double itkNotUsed(fc),
double fa, double functionValueOfb, double fc,
double * extX, double * extVal)
{
PowellOptimizer::ParametersType tempCoord( m_SpaceDimension );
this->BracketedLineOptimize( ax, bx, cx, fa, functionValueOfb, fc, extX, extVal, tempCoord);
}

void
PowellOptimizer
::BracketedLineOptimize(double ax, double bx, double cx,
double itkNotUsed(fa), double functionValueOfb, double itkNotUsed(fc),
double * extX, double * extVal,
ParametersType & tempCoord)
{
double x;
double v = 0.0;
Expand Down Expand Up @@ -389,7 +415,7 @@ ::BracketedLineOptimize(double ax, double bx, double cx,

double functionValueOft;

functionValueOft = this->GetLineValue(t);
functionValueOft = this->GetLineValue(t, tempCoord);

if( functionValueOft <= functionValueOfX )
{
Expand Down Expand Up @@ -459,13 +485,16 @@ ::StartOptimization()
m_SpaceDimension = m_CostFunction->GetNumberOfParameters();
m_LineOrigin.set_size(m_SpaceDimension);
m_LineDirection.set_size(m_SpaceDimension);
this->m_CurrentPosition.set_size(m_SpaceDimension);

vnl_matrix<double> xi(m_SpaceDimension, m_SpaceDimension);
vnl_vector<double> xit(m_SpaceDimension);
xi.set_identity();
xit.fill(0);
xit[0] = 1;

PowellOptimizer::ParametersType tempCoord(m_SpaceDimension);

PowellOptimizer::ParametersType p(m_SpaceDimension);
PowellOptimizer::ParametersType pt(m_SpaceDimension);
PowellOptimizer::ParametersType ptt(m_SpaceDimension);
Expand All @@ -479,7 +508,7 @@ ::StartOptimization()

xx = 0;
this->SetLine(p, xit);
fx = this->GetLineValue(0);
fx = this->GetLineValue(0, tempCoord);

for (m_CurrentIteration = 0;
m_CurrentIteration <= m_MaximumIteration;
Expand All @@ -502,8 +531,8 @@ ::StartOptimization()
ax = 0.0;
fa = fx;
xx = m_StepLength;
this->LineBracket( &ax, &xx, &bx, &fa, &fx, &fb);
this->BracketedLineOptimize(ax, xx, bx, fa, fx, fb, &xx, &fx);
this->LineBracket( &ax, &xx, &bx, &fa, &fx, &fb, tempCoord);
this->BracketedLineOptimize(ax, xx, bx, fa, fx, fb, &xx, &fx, tempCoord);
this->SetCurrentLinePoint(xx, fx);
p = this->GetCurrentPosition();

Expand All @@ -529,7 +558,7 @@ ::StartOptimization()
}

this->SetLine(ptt, xit);
fptt = this->GetLineValue(0);
fptt = this->GetLineValue(0, tempCoord);
if (fptt < fp)
{
double t = 2.0 * (fp - 2.0*fx + fptt)
Expand All @@ -542,8 +571,8 @@ ::StartOptimization()
ax = 0.0;
fa = fx;
xx = 1;
this->LineBracket( &ax, &xx, &bx, &fa, &fx, &fb);
this->BracketedLineOptimize(ax, xx, bx, fa, fx, fb, &xx, &fx);
this->LineBracket( &ax, &xx, &bx, &fa, &fx, &fb, tempCoord);
this->BracketedLineOptimize(ax, xx, bx, fa, fx, fb, &xx, &fx, tempCoord);
this->SetCurrentLinePoint(xx, fx);
p = this->GetCurrentPosition();

Expand Down
20 changes: 19 additions & 1 deletion Code/Numerics/itkPowellOptimizer.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -145,6 +145,8 @@ class ITK_EXPORT PowellOptimizer:
* Line origin and distances are set via SetLine */
double GetLineValue(double x) const;

double GetLineValue(double x, ParametersType & tempCoord) const;

/** Set the given scalar step distance (x) and function value (fx) as the
* "best-so-far" optimizer values. */
void SetCurrentLinePoint(double x, double fx);
Expand All @@ -169,6 +171,10 @@ class ITK_EXPORT PowellOptimizer:
virtual void LineBracket(double *ax, double *bx, double *cx,
double *fa, double *fb, double *fc);

virtual void LineBracket(double *ax, double *bx, double *cx,
double *fa, double *fb, double *fc,
ParametersType & tempCoord);

/** Given a bracketing triple of points and their function values, returns
* a bounded extreme. These values are in parameter space, along the
* current line and wrt the current origin set via SetLine. Optimization
Expand All @@ -178,8 +184,20 @@ class ITK_EXPORT PowellOptimizer:
double fa, double fb, double fc,
double * extX, double * extVal);

virtual void BracketedLineOptimize(double ax, double bx, double cx,
double fa, double fb, double fc,
double * extX, double * extVal,
ParametersType & tempCoord);

itkGetMacro(SpaceDimension, unsigned int);
itkSetMacro(SpaceDimension, unsigned int);
void SetSpaceDimension( unsigned int dim )
{
this->m_SpaceDimension = dim;
this->m_LineDirection.set_size( dim );
this->m_LineOrigin.set_size( dim );
this->m_CurrentPosition.set_size( dim );
this->Modified();
}

itkSetMacro(CurrentIteration, unsigned int);

Expand Down

0 comments on commit a7b0a73

Please sign in to comment.