Commit 28e74220 authored by Roman Tolchenov's avatar Roman Tolchenov
Browse files

Re #5893. Added option to calculate derivatives numerically

 to composite function. Using parameter scaling in LevenbergMarquardtMD.
parent 7b5880f5
......@@ -54,7 +54,7 @@ class MANTID_API_DLL CompositeFunction : public virtual IFunction
{
public:
/// Default constructor
CompositeFunction(): m_nParams(0){}
CompositeFunction(): m_nParams(0),m_useNumericDerivatives(false){}
/// Copy contructor
//CompositeFunction(const CompositeFunction&);
///Assignment operator
......@@ -161,14 +161,14 @@ public:
void replaceFunctionPtr(const IFunction_sptr f_old,IFunction_sptr f_new);
/// Get the function index
std::size_t functionIndex(std::size_t i)const;
/// Get the function index
//std::size_t functionIndexActive(std::size_t i)const;
/// Returns the index of parameter i as it declared in its function
size_t parameterLocalIndex(size_t i)const;
/// Returns the name of parameter i as it declared in its function
std::string parameterLocalName(size_t i)const;
/// Check the function.
void checkFunction();
/// Enable/disable numeric derivative calculation
void useNumericDerivatives( bool yes ) const;
/// Returns the number of attributes associated with the function
virtual size_t nLocalAttributes()const {return 0;}
......@@ -217,6 +217,8 @@ private:
size_t m_nParams;
/// Function counter to be used in nextConstraint
mutable size_t m_iConstraintFunction;
/// Flag set to use numerical derivatives
mutable bool m_useNumericDerivatives;
};
///shared pointer to the composite function base class
......
......@@ -121,10 +121,17 @@ void CompositeFunction::function(const FunctionDomain& domain, FunctionValues& v
*/
void CompositeFunction::functionDeriv(const FunctionDomain& domain, Jacobian& jacobian)
{
for(size_t iFun = 0; iFun < nFunctions(); ++iFun)
if ( m_useNumericDerivatives )
{
calNumericalDeriv(domain, jacobian);
}
else
{
PartialJacobian J(&jacobian,paramOffset(iFun));
getFunction(iFun)->functionDeriv(domain,J);
for(size_t iFun = 0; iFun < nFunctions(); ++iFun)
{
PartialJacobian J(&jacobian,paramOffset(iFun));
getFunction(iFun)->functionDeriv(domain,J);
}
}
}
......@@ -601,18 +608,6 @@ size_t CompositeFunction::functionIndex(std::size_t i)const
return m_IFunction[i];
}
/**
* Get the index of the function to which parameter i belongs
* @param i :: The active parameter index
* @return active function index of the requested parameter
*/
//size_t CompositeFunction::functionIndexActive(std::size_t i)const
//{
// if( i >= nParams() )
// throw std::out_of_range("Function parameter index out of range.");
// return m_IFunctionActive[i];
//}
/**
* @param varName :: The variable name which may contain function index ( [f<index.>]name )
* @param index :: Receives function index or throws std::invalid_argument
......@@ -820,5 +815,15 @@ IFunction_sptr CompositeFunction::getContainingFunction(const ParameterReference
return IFunction_sptr();
}
/**
* Enable/disable numeric derivative calculation.
* @param yes :: Set to true to use numeric derivative calculation.
*/
void CompositeFunction::useNumericDerivatives( bool yes ) const
{
m_useNumericDerivatives = yes;
}
} // namespace API
} // namespace Mantid
......@@ -28,9 +28,16 @@ FuncMinimizerFactoryImpl::FuncMinimizerFactoryImpl() : Kernel::DynamicFactory<IF
*/
boost::shared_ptr<IFuncMinimizer> FuncMinimizerFactoryImpl::createMinimizer(const std::string& str) const
{
// check if there are any properties defined - look for a comma
if ( str.find(',') == std::string::npos )
{// no properties - create minimizer and return
return create( str );
}
// parse the string
Expression parser;
parser.parse( str );
parser.toList(); // make it a list even if it has 1 item
parser.toList(); // make sure it is a list
const size_t n = parser.size();
if ( n == 0 )
{
......@@ -43,7 +50,7 @@ boost::shared_ptr<IFuncMinimizer> FuncMinimizerFactoryImpl::createMinimizer(cons
const std::string type = parser[0].str();
auto minimizer = create( type );
// set the properties if there are any
// set the properties
for(size_t i = 1; i < n; ++i)
{
auto& param = parser[i];
......
......@@ -1386,6 +1386,7 @@ bool LeBailFit::fitLeBailFunction(size_t workspaceindex, std::map<std::string, P
g_log.debug() << "[Before Fit] Function To Fit: " << mLeBailFunction->asString() << std::endl;
// b) Set property
mLeBailFunction->useNumericDerivatives( true );
fitalg->setProperty("Function", boost::shared_ptr<API::IFunction>(mLeBailFunction));
fitalg->setProperty("InputWorkspace", dataWS);
fitalg->setProperty("WorkspaceIndex", int(workspaceindex));
......
......@@ -25,7 +25,7 @@ namespace CurveFitting
DECLARE_FUNCMINIMIZER(LevenbergMarquardtMDMinimizer,Levenberg-MarquardtMD)
// Get a reference to the logger
Kernel::Logger& LevenbergMarquardtMDMinimizer::g_log = Kernel::Logger::get("LevenbergMarquardtMDMinimizer");
Kernel::Logger& LevenbergMarquardtMDMinimizer::g_log = Kernel::Logger::get("Lev-MarqMD");
/// Constructor
LevenbergMarquardtMDMinimizer::LevenbergMarquardtMDMinimizer():
......@@ -36,6 +36,8 @@ m_mu(0),
m_nu(2.0),
m_rho(1.0)
{
declareProperty("MuMax",1e6,"Maximum value of mu.");
declareProperty("Debug",false,"Turn on the debug output.");
}
/// Initialize minimizer, i.e. pass a function to minimize.
......@@ -54,7 +56,8 @@ void LevenbergMarquardtMDMinimizer::initialize(API::ICostFunction_sptr function)
/// Do one iteration.
bool LevenbergMarquardtMDMinimizer::iterate()
{
const bool debug = false;
const bool debug = getProperty("Debug");
const double muMax = getProperty("MuMax");
if ( !m_leastSquares )
{
......@@ -64,7 +67,13 @@ bool LevenbergMarquardtMDMinimizer::iterate()
if ( n == 0 )
{
m_errorString = "No parameters to fit";
m_errorString = "No parameters to fit.";
return false;
}
if ( m_mu > muMax )
{
m_errorString = "Failed to converge.";
return false;
}
......@@ -88,7 +97,7 @@ bool LevenbergMarquardtMDMinimizer::iterate()
if (debug)
{
std::cerr << "mu=" << m_mu << std::endl;
g_log.warning() << "mu=" << m_mu << std::endl;
}
if (m_D.empty())
......@@ -97,8 +106,11 @@ bool LevenbergMarquardtMDMinimizer::iterate()
}
// copy the hessian
GSLMatrix H(m_leastSquares->getHessian()/*m_hessian*/);
GSLVector dd(m_leastSquares->getDeriv()/*m_der*/);
GSLMatrix H( m_leastSquares->getHessian() );
GSLVector dd( m_leastSquares->getDeriv() );
// scaling factors
std::vector<double> sf(n);
for(size_t i = 0; i < n; ++i)
{
......@@ -107,17 +119,41 @@ bool LevenbergMarquardtMDMinimizer::iterate()
m_D[i] = d;
double tmp = H.get(i,i) + m_mu * d;
H.set(i,i,tmp);
sf[i] = sqrt( tmp );
if ( tmp == 0.0 )
{
m_errorString = "Singular matrix.";
return false;
}
}
// apply scaling
for(size_t i = 0; i < n; ++i)
{
double d = dd.get(i);
dd.set( i, d / sf[i] );
for(size_t j = i; j < n; ++j)
{
const double f = sf[i] * sf[j];
double tmp = H.get(i,j);
H.set(i,j,tmp/f);
if ( i != j )
{
tmp = H.get(j,i);
H.set(j,i,tmp/f);
}
}
}
if (debug && m_rho > 0)
{
std::cerr << "H:\n" << H ;
std::cerr << "-----------------------------\n";
for(size_t j = 0; j < n; ++j) {std::cerr << dd.get(j) << ' '; } std::cerr << std::endl;
std::cerr << "det=" << H.det() << std::endl;
g_log.warning() << "H:\n" << H ;
g_log.warning() << "-----------------------------\n";
for(size_t j = 0; j < n; ++j) {g_log.warning() << dd.get(j) << ' '; } g_log.warning() << std::endl;
g_log.warning() << "det=" << H.det() << std::endl;
}
/// Parameter corrections
// Parameter corrections
GSLVector dx(n);
// To find dx solve the system of linear equations H * dx == -m_der
dd *= -1.0;
......@@ -125,7 +161,19 @@ bool LevenbergMarquardtMDMinimizer::iterate()
if (debug)
{
for(size_t j = 0; j < n; ++j) {std::cerr << dx.get(j) << ' '; } std::cerr << std::endl << std::endl;
g_log.warning() << "\nScaling:" << std::endl;
for(size_t j = 0; j < n; ++j) {g_log.warning() << sf[j] << ' '; } g_log.warning() << std::endl;
g_log.warning() << "Corrections:" << std::endl;
for(size_t j = 0; j < n; ++j) {g_log.warning() << dx.get(j) << ' '; } g_log.warning() << std::endl << std::endl;
}
// restore scaling
for(size_t i = 0; i < n; ++i)
{
double d = dx.get(i);
dx.set( i, d / sf[i] );
d = dd.get(i);
dd.set( i, d * sf[i] );
}
// save previous state
......@@ -137,7 +185,7 @@ bool LevenbergMarquardtMDMinimizer::iterate()
m_leastSquares->setParameter(i,d);
if (debug)
{
std::cerr << "P" << i << ' ' << d << std::endl;
g_log.warning() << "P" << i << ' ' << d << std::endl;
}
}
m_leastSquares->getFittingFunction()->applyTies();
......@@ -155,8 +203,8 @@ bool LevenbergMarquardtMDMinimizer::iterate()
if (debug)
{
static size_t iter = 0;
std::cerr << "iter " << ++iter << std::endl;
std::cerr << "F " << m_F << ' ' << F1 << ' ' << dL << std::endl;
g_log.warning() << "iter " << ++iter << std::endl;
g_log.warning() << "F " << m_F << ' ' << F1 << ' ' << dL << std::endl;
}
// Try the stop condition
......@@ -196,7 +244,7 @@ bool LevenbergMarquardtMDMinimizer::iterate()
}
if (debug)
{
std::cerr << "rho=" << m_rho << std::endl;
g_log.warning() << "rho=" << m_rho << std::endl;
}
if (m_rho > 0)
......@@ -212,7 +260,7 @@ bool LevenbergMarquardtMDMinimizer::iterate()
m_nu = 2.0;
m_F = F1;
if (debug)
std::cerr << "times " << m_rho << std::endl;
g_log.warning() << "times " << m_rho << std::endl;
// drop saved state, accept new parameters
m_leastSquares->drop();
}
......@@ -225,6 +273,7 @@ bool LevenbergMarquardtMDMinimizer::iterate()
m_F = m_leastSquares->val();
}
//m_leastSquares->drop();
return true;
}
......
......@@ -1042,7 +1042,19 @@ std::string FitPropertyBrowser::minimizer(bool withProperties)const
{
foreach(QtProperty* prop,m_minimizerProperties)
{
minimStr += "," + prop->propertyName() + "=" + QString::number( m_doubleManager->value(prop) );
minimStr += "," + prop->propertyName() + "=";
if ( prop->propertyManager() == m_doubleManager )
{
minimStr += QString::number( m_doubleManager->value(prop) );
}
else if ( prop->propertyManager() == m_boolManager )
{
minimStr += QString::number( m_boolManager->value(prop) );
}
else
{
minimStr += m_stringManager->value(prop);
}
}
}
return minimStr.toStdString();
......@@ -2952,7 +2964,24 @@ void FitPropertyBrowser::minimizerChanged()
for(auto it = properties.begin(); it != properties.end(); ++it)
{
QString propName = QString::fromStdString( (**it).name() );
QtProperty* prop = this->addDoubleProperty( propName );
QtProperty* prop = NULL;
if ( auto prp = dynamic_cast<Mantid::Kernel::PropertyWithValue<bool>* >(*it) )
{
prop = m_boolManager->addProperty( propName );
bool val = *prp;
m_boolManager->setValue( prop, val );
}
else if ( auto prp = dynamic_cast<Mantid::Kernel::PropertyWithValue<double>* >(*it) )
{
prop = this->addDoubleProperty( propName );
double val = *prp;
m_doubleManager->setValue( prop, val );
}
else
{
prop = m_stringManager->addProperty( propName );
QString val = QString::fromStdString( prp->value() );
}
m_settingsGroup->property()->addSubProperty( prop );
m_minimizerProperties.append( prop );
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment