From 0cd3ae9e217d9d1736083b2b4e1bbed4cac0b2b3 Mon Sep 17 00:00:00 2001 From: Vickie Lynch Date: Tue, 11 Nov 2014 08:47:09 -0500 Subject: [PATCH] Refs #10199 use same cost function as ISAW --- .../Crystal/src/OptimizeCrystalPlacement.cpp | 12 + .../Framework/Crystal/src/PeakHKLErrors.cpp | 47 +- .../Framework/CurveFitting/CMakeLists.txt | 2 + .../inc/MantidCurveFitting/CostFuncSimplex.h | 115 ++++ .../inc/MantidCurveFitting/SeqDomain.h | 3 + .../CurveFitting/src/CostFuncSimplex.cpp | 539 ++++++++++++++++++ .../Framework/CurveFitting/src/SeqDomain.cpp | 45 +- .../MantidGeometry/Crystal/IndexingUtils.h | 4 +- .../Geometry/src/Crystal/IndexingUtils.cpp | 4 +- 9 files changed, 746 insertions(+), 25 deletions(-) create mode 100644 Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/CostFuncSimplex.h create mode 100644 Code/Mantid/Framework/CurveFitting/src/CostFuncSimplex.cpp diff --git a/Code/Mantid/Framework/Crystal/src/OptimizeCrystalPlacement.cpp b/Code/Mantid/Framework/Crystal/src/OptimizeCrystalPlacement.cpp index 793c6a5d9012..cad2e925991b 100644 --- a/Code/Mantid/Framework/Crystal/src/OptimizeCrystalPlacement.cpp +++ b/Code/Mantid/Framework/Crystal/src/OptimizeCrystalPlacement.cpp @@ -355,6 +355,10 @@ namespace Mantid fit_alg->setProperty( "MaxIterations" , 60 ); + fit_alg->setProperty( "Minimizer" , "Simplex" ); + + fit_alg->setProperty( "CostFunction" , "Simplex" ); + fit_alg->setProperty( "Constraints" , Constr ); fit_alg->setProperty( "InputWorkspace" , mwkspc ); @@ -380,6 +384,14 @@ namespace Mantid Ties +="GonRotx=0.0,GonRoty=0.0,GonRotz=0.0"; } + else + { + if( !Ties.empty()) + Ties +=","; + + Ties +="GonRotx=0.0,GonRotz=0.0"; + + } if( !Ties.empty()) fit_alg->setProperty("Ties",Ties); diff --git a/Code/Mantid/Framework/Crystal/src/PeakHKLErrors.cpp b/Code/Mantid/Framework/Crystal/src/PeakHKLErrors.cpp index 6a3eb3f0ffa3..2b8f27292aa8 100644 --- a/Code/Mantid/Framework/Crystal/src/PeakHKLErrors.cpp +++ b/Code/Mantid/Framework/Crystal/src/PeakHKLErrors.cpp @@ -21,12 +21,14 @@ #include "MantidKernel/Matrix.h" #include "MantidKernel/Quat.h" #include "MantidKernel/V3D.h" +#include "MantidGeometry/Crystal/IndexingUtils.h" #include #include using namespace Mantid::DataObjects; using namespace Mantid::API; using namespace Mantid::Kernel; +using namespace Mantid::Geometry; using Mantid::Geometry::CompAssembly; using Mantid::Geometry::IObjComponent_const_sptr; using Mantid::Geometry::IComponent_const_sptr; @@ -364,11 +366,7 @@ namespace Mantid std::map > RunNum2GonMatrixMap; getRun2MatMap( Peaks, OptRuns, RunNum2GonMatrixMap); - const DblMatrix & UBx = Peaks->sample().getOrientedLattice().getUB(); - - DblMatrix UBinv( UBx ); - UBinv.Invert(); - UBinv /= ( 2 * M_PI ); + DblMatrix UB = Peaks->sample().getOrientedLattice().getUB(); double GonRotx =getParameter("GonRotx"); double GonRoty =getParameter("GonRoty"); @@ -379,6 +377,7 @@ namespace Mantid double ChiSqTot = 0.0; + std::vector q_vectors; for( size_t i = 0; i< nData; i+= 3 ) { int peakNum = (int)( .5 + xValues[i] ); @@ -400,27 +399,35 @@ namespace Mantid peak.setGoniometerMatrix( GonRot*peak.getGoniometerMatrix()); } - V3D hkl = UBinv * peak.getQSampleFrame(); - - - - for( int k = 0; k<3; k++ ) - { - double d1 = hkl[k] - floor( hkl[k] ); - if( d1>.5 ) d1 = d1 - 1; - if( d1 < -.5 ) d1 = d1 + 1; - - out[i+k] = d1; - ChiSqTot += d1*d1; - } - + q_vectors.push_back( peak.getQSampleFrame()); } + std::vector sigabc(1);// only need chisq + std::vector miller_ind; + std::vector indexed_qs; + double fit_error; + miller_ind.reserve( q_vectors.size() ); + indexed_qs.reserve( q_vectors.size() ); + double tolerance = 0.12; + IndexingUtils::GetIndexedPeaks( UB, q_vectors, tolerance, + miller_ind, indexed_qs, fit_error ); + ChiSqTot = IndexingUtils::Optimize_UB(UB, miller_ind,indexed_qs,out); + + int num_indexed = IndexingUtils::NumberIndexed(UB, q_vectors, tolerance); + /*for( size_t p = 0 ; p < nData; p++ ) + { + out[p] *= 0.1*static_cast(num_indexed)*static_cast(nData); + out[p] -= static_cast(num_indexed)/static_cast(nData); + }*/ + g_log.notice() << "New UB will index " << num_indexed << " Peaks out of " << nData + << " with tolerance of " << std::setprecision(3) << std::setw(5) << tolerance + << "\n"; g_log.debug() << "------------------------Function-----------------------------------------------"<. + + File change history is stored at: . + Code Documentation is available at: +*/ +class DLLExport CostFuncSimplex : public CostFuncFitting +{ +public: + /// Constructor + CostFuncSimplex(); + /// Virtual destructor + virtual ~CostFuncSimplex() {} + + /// Get name of minimizer + virtual std::string name() const { return "Simplex";} + + /// Get short name of minimizer - useful for say labels in guis + virtual std::string shortName() const {return "Chi";}; + + /// Calculate value of cost function + virtual double val() const; + + /// Calculate the derivatives of the cost function + /// @param der :: Container to output the derivatives + virtual void deriv(std::vector& der) const; + + /// Calculate the value and the derivatives of the cost function + /// @param der :: Container to output the derivatives + /// @return :: The value of the function + virtual double valAndDeriv(std::vector& der) const; + + virtual double valDerivHessian(bool evalFunction = true, bool evalDeriv = true, bool evalHessian = true) const; + const GSLVector& getDeriv() const; + const GSLMatrix& getHessian() const; + void push(); + void pop(); + void drop(); + + void setParameters(const GSLVector& params); + void getParameters(GSLVector& params) const; +protected: + + virtual void calActiveCovarianceMatrix(GSLMatrix& covar, double epsrel = 1e-8); + + void addVal( + API::FunctionDomain_sptr domain, + API::FunctionValues_sptr values + )const; + void addValDerivHessian( + API::IFunction_sptr function, + API::FunctionDomain_sptr domain, + API::FunctionValues_sptr values, + bool evalFunction = true, bool evalDeriv = true, bool evalHessian = true) const; + + /// Get mapped weights from FunctionValues + virtual std::vector getFitWeights(API::FunctionValues_sptr values) const; + + /// Flag to include constraint in cost function value + bool m_includePenalty; + + mutable double m_value; + mutable GSLVector m_der; + mutable GSLMatrix m_hessian; + + mutable bool m_pushed; + mutable double m_pushedValue; + mutable GSLVector m_pushedParams; + + friend class SeqDomain; + friend class ParDomain; + + double m_factor; +}; + +} // namespace CurveFitting +} // namespace Mantid + +#endif /*MANTID_CURVEFITTING_COSTFUNCSIMPLEX_H_*/ diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/SeqDomain.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/SeqDomain.h index e5d1fb776e12..44e451c788aa 100644 --- a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/SeqDomain.h +++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/SeqDomain.h @@ -10,6 +10,7 @@ #include "MantidAPI/IDomainCreator.h" #include "MantidCurveFitting/CostFuncLeastSquares.h" #include "MantidCurveFitting/CostFuncRwp.h" +#include "MantidCurveFitting/CostFuncSimplex.h" #include #include @@ -65,6 +66,8 @@ class MANTID_CURVEFITTING_DLL SeqDomain: public API::FunctionDomain void rwpVal(const CostFuncRwp& rwp); /// Calculate the value, first and second derivatives of a RWP cost function void rwpValDerivHessian(const CostFuncRwp& rwp, bool evalFunction, bool evalDeriv, bool evalHessian); + virtual void simplexVal(const CostFuncSimplex& simplex); + virtual void simplexValDerivHessian(const CostFuncSimplex& leastSquares, bool evalFunction, bool evalDeriv, bool evalHessian); /// Create an instance of SeqDomain in one of two forms: either SeqDomain for sequential domain creation /// or ParDomain for parallel calculations diff --git a/Code/Mantid/Framework/CurveFitting/src/CostFuncSimplex.cpp b/Code/Mantid/Framework/CurveFitting/src/CostFuncSimplex.cpp new file mode 100644 index 000000000000..308de03f4158 --- /dev/null +++ b/Code/Mantid/Framework/CurveFitting/src/CostFuncSimplex.cpp @@ -0,0 +1,539 @@ +//---------------------------------------------------------------------- +// Includes +//---------------------------------------------------------------------- +#include "MantidCurveFitting/CostFuncSimplex.h" +#include "MantidCurveFitting/Jacobian.h" +#include "MantidCurveFitting/SeqDomain.h" +#include "MantidAPI/IConstraint.h" +#include "MantidAPI/CompositeDomain.h" +#include "MantidAPI/FunctionValues.h" +#include "MantidKernel/Logger.h" +#include "MantidKernel/MultiThreaded.h" + +#include + +namespace Mantid +{ +namespace CurveFitting +{ + namespace + { + /// static logger + Kernel::Logger g_log("CostFuncSimplex"); + } + + +DECLARE_COSTFUNCTION(CostFuncSimplex,Simplex) + +/** + * Constructor + */ +CostFuncSimplex::CostFuncSimplex() : CostFuncFitting(), + m_includePenalty(true), + m_value(0), + m_pushed(false), + m_factor(0.5) + { + } + +/** Calculate value of cost function + * @return :: The value of the function + */ +double CostFuncSimplex::val() const +{ + if ( !m_dirtyVal ) return m_value; + + checkValidity(); + + m_value = 0.0; + + auto seqDomain = boost::dynamic_pointer_cast(m_domain); + + if (seqDomain) + { + seqDomain->simplexVal(*this); + } + else + { + if (!m_values) + { + throw std::runtime_error("Simplex: undefined FunctionValues."); + } + addVal(m_domain,m_values); + } + + // add penalty + if (m_includePenalty) + { + for(size_t i=0;inParams();++i) + { + if ( !m_function->isActive(i) ) continue; + API::IConstraint* c = m_function->getConstraint(i); + if (c) + { + m_value += c->check(); + } + } + } + + m_dirtyVal = false; + return m_value; +} + +/** + * Add a contribution to the cost function value from the fitting function evaluated on a particular domain. + * @param domain :: A domain + * @param values :: Values + */ +void CostFuncSimplex::addVal(API::FunctionDomain_sptr domain, API::FunctionValues_sptr values)const +{ + m_function->function(*domain,*values); + size_t ny = values->size(); + + double retVal = 0.0; + + std::vector weights = getFitWeights(values); + + for (size_t i = 0; i < ny; i++) + { + // double val = ( values->getCalculated(i) - values->getFitData(i) ) * values->getFitWeight(i); + double val = ( values->getCalculated(i) - values->getFitData(i) ) * weights[i]; + retVal += val; + g_log.debug() << values->getCalculated(i)<<" "<getFitData(i)<<" "<& der) const +{ + valDerivHessian(false,true,false); + + if (der.size() != nParams()) + { + der.resize(nParams()); + } + for(size_t i = 0; i < nParams(); ++i) + { + der[i] = m_der.get(i); + } +} + +/** Calculate the value and the derivatives of the cost function + * @param der :: Container to output the derivatives + * @return :: The value of the function + */ +double CostFuncSimplex::valAndDeriv(std::vector& der) const +{ + valDerivHessian(true,true,false); + + if (der.size() != nParams()) + { + der.resize(nParams()); + } + for(size_t i = 0; i < nParams(); ++i) + { + der[i] = m_der.get(i); + } + return m_value; +} + +/** Calculate the value and the first and second derivatives of the cost function + * @param evalFunction :: If false cost function isn't evaluated and returned value (0.0) should be ignored. + * It is for efficiency reasons. + * @param evalDeriv :: flag for evaluation of the first derivatives + * @param evalHessian :: flag for evaluation of the second derivatives + */ +double CostFuncSimplex::valDerivHessian(bool evalFunction, bool evalDeriv, bool evalHessian) const +{ + if (m_pushed || !evalDeriv) + { + return val(); + } + + if (!m_dirtyVal && !m_dirtyDeriv && !m_dirtyHessian) return m_value; + if (m_dirtyVal) evalFunction = true; + + checkValidity(); + + if (evalFunction) + { + m_value = 0.0; + } + if (evalDeriv) + { + m_der.resize(nParams()); + m_der.zero(); + } + if (evalHessian) + { + m_hessian.resize(nParams(),nParams()); + m_hessian.zero(); + } + + auto seqDomain = boost::dynamic_pointer_cast(m_domain); + + if (seqDomain) + { + seqDomain->simplexValDerivHessian(*this,evalFunction,evalDeriv,evalHessian); + } + else + { + if (!m_values) + { + throw std::runtime_error("Simplex: undefined FunctionValues."); + } + addValDerivHessian(m_function,m_domain,m_values,evalFunction,evalDeriv,evalHessian); + } + + // Add constraints penalty + size_t np = m_function->nParams(); + if (evalFunction) + { + if (m_includePenalty) + { + for(size_t i = 0; i < np; ++i) + { + API::IConstraint* c = m_function->getConstraint(i); + if (c) + { + m_value += c->check(); + } + } + } + m_dirtyVal = false; + } + + if (evalDeriv) + { + if (m_includePenalty) + { + size_t i = 0; + for(size_t ip = 0; ip < np; ++ip) + { + if ( !m_function->isActive(ip) ) continue; + API::IConstraint* c = m_function->getConstraint(ip); + if (c) + { + double d = m_der.get(i) + c->checkDeriv(); + m_der.set(i,d); + } + ++i; + } + } + m_dirtyDeriv = false; + } + + if (evalDeriv) + { + if (m_includePenalty) + { + size_t i = 0; + for(size_t ip = 0; ip < np; ++ip) + { + if ( !m_function->isActive(ip) ) continue; + API::IConstraint* c = m_function->getConstraint(ip); + if (c) + { + double d = m_hessian.get(i,i) + c->checkDeriv2(); + m_hessian.set(i,i,d); + } + ++i; + } + } + // clear the dirty flag if hessian was actually calculated + m_dirtyHessian = m_hessian.isEmpty(); + } + + return m_value; +} + +/** + * Update the cost function, derivatives and hessian by adding values calculated + * on a domain. + * @param function :: Function to use to calculate the value and the derivatives + * @param domain :: The domain. + * @param values :: The fit function values + * @param evalFunction :: Flag to evaluate the function + * @param evalDeriv :: Flag to evaluate the derivatives + * @param evalHessian :: Flag to evaluate the Hessian + */ +void CostFuncSimplex::addValDerivHessian( + API::IFunction_sptr function, + API::FunctionDomain_sptr domain, + API::FunctionValues_sptr values, + bool evalFunction , bool evalDeriv, bool evalHessian) const +{ + UNUSED_ARG(evalDeriv); + size_t np = function->nParams(); // number of parameters + size_t ny = domain->size(); // number of data points + Jacobian jacobian(ny,np); + if (evalFunction) + { + function->function(*domain,*values); + } + function->functionDeriv(*domain,jacobian); + + size_t iActiveP = 0; + double fVal = 0.0; + if (g_log.is(Kernel::Logger::Priority::PRIO_DEBUG)) + { + g_log.debug() << "Jacobian:\n"; + for(size_t i = 0; i < ny; ++i) + { + for(size_t ip = 0; ip < np; ++ip) + { + if ( !m_function->isActive(ip) ) continue; + g_log.debug() << jacobian.get(i,ip) << ' '; + } + g_log.debug() << "\n"; + } + } + + std::vector weights = getFitWeights(values); + + for(size_t ip = 0; ip < np; ++ip) + { + if ( !function->isActive(ip) ) continue; + double d = 0.0; + for(size_t i = 0; i < ny; ++i) + { + double calc = values->getCalculated(i); + double obs = values->getFitData(i); + double w = weights[i]; + double y = ( calc - obs ) * w; + d += y * jacobian.get(i,ip) * w; + if (iActiveP == 0 && evalFunction) + { + fVal += y * y; + } + } + PARALLEL_CRITICAL(der_set) + { + double der = m_der.get(iActiveP); + m_der.set(iActiveP, der + d); + } + //std::cerr << "der " << ip << ' ' << der[iActiveP] << std::endl; + ++iActiveP; + } + + if (evalFunction) + { + PARALLEL_ATOMIC + m_value += 0.5 * fVal; + } + + if (!evalHessian) return; + + //size_t na = m_der.size(); // number of active parameters + + size_t i1 = 0; // active parameter index + for(size_t i = 0; i < np; ++i) // over parameters + { + if ( !function->isActive(i) ) continue; + size_t i2 = 0; // active parameter index + for(size_t j = 0; j <= i; ++j) // over ~ half of parameters + { + if ( !function->isActive(j) ) continue; + double d = 0.0; + for(size_t k = 0; k < ny; ++k) // over fitting data + { + // double w = values->getFitWeight(k); + double w = weights[k]; + d += jacobian.get(k,i) * jacobian.get(k,j) * w * w; + } + PARALLEL_CRITICAL(hessian_set) + { + double h = m_hessian.get(i1,i2); + m_hessian.set(i1,i2, h + d); + //std::cerr << "hess " << i1 << ' ' << i2 << std::endl; + if (i1 != i2) + { + m_hessian.set(i2,i1,h + d); + } + } + ++i2; + } + ++i1; + } +} + +std::vector CostFuncSimplex::getFitWeights(API::FunctionValues_sptr values) const +{ + std::vector weights(values->size()); + for(size_t i = 0; i < weights.size(); ++i) { + weights[i] = values->getFitWeight(i); + } + + return weights; +} + +/** + * Return cached or calculate the drivatives. + */ +const GSLVector& CostFuncSimplex::getDeriv() const +{ + if (m_pushed) + { + return m_der; + } + if (m_dirtyVal || m_dirtyDeriv || m_dirtyHessian) + { + valDerivHessian(); + } + return m_der; +} + +/** + * Return cached or calculate the Hessian. + */ +const GSLMatrix& CostFuncSimplex::getHessian() const +{ + if (m_pushed) + { + return m_hessian; + } + if (m_dirtyVal || m_dirtyDeriv || m_dirtyHessian) + { + valDerivHessian(); + } + return m_hessian; +} + +/** + * Save current parameters, derivatives and hessian. + */ +void CostFuncSimplex::push() +{ + if (m_pushed) + { + throw std::runtime_error("Simplex: double push."); + } + // make sure we are not dirty + m_pushedValue = valDerivHessian(); + getParameters(m_pushedParams); + m_pushed = true; +} + +/** + * Restore saved parameters, derivatives and hessian. + */ +void CostFuncSimplex::pop() +{ + if ( !m_pushed ) + { + throw std::runtime_error("Simplex: empty stack."); + } + setParameters(m_pushedParams); + m_value = m_pushedValue; + m_pushed = false; + m_dirtyVal = false; + m_dirtyDeriv = false; + m_dirtyHessian = false; +} + +/** + * Discard saved parameters, derivatives and hessian. + */ +void CostFuncSimplex::drop() +{ + if ( !m_pushed ) + { + throw std::runtime_error("Simplex: empty stack."); + } + m_pushed = false; + setDirty(); +} + +/** + * Copy the parameter values from a GSLVector. + * @param params :: A vector to copy the parameters from + */ +void CostFuncSimplex::setParameters(const GSLVector& params) +{ + if (nParams() != params.size()) + { + throw std::runtime_error("Parameter vector has wrong size in CostFuncSimplex."); + } + for(size_t i = 0; i < nParams(); ++i) + { + setParameter(i,params.get(i)); + } + m_function->applyTies(); +} + +/** + * Copy the parameter values to a GSLVector. + * @param params :: A vector to copy the parameters to + */ +void CostFuncSimplex::getParameters(GSLVector& params) const +{ + if (params.size() != nParams()) + { + params.resize(nParams()); + } + for(size_t i = 0; i < nParams(); ++i) + { + params.set(i,getParameter(i)); + } +} + +/** + * Calculates covariance matrix for fitting function's active parameters. + * @param covar :: Output cavariance matrix. + * @param epsrel :: Tolerance. + */ +void CostFuncSimplex::calActiveCovarianceMatrix(GSLMatrix& covar, double epsrel) +{ + UNUSED_ARG(epsrel); + if (m_hessian.isEmpty()) + { + valDerivHessian(); + } + if(g_log.is(Kernel::Logger::Priority::PRIO_INFORMATION)) + { + g_log.information() << "== Hessian (H) ==\n"; + std::ios::fmtflags prevState = g_log.information().flags(); + g_log.information() << std::left << std::fixed; + for(size_t i = 0; i < m_hessian.size1(); ++i) + { + for(size_t j = 0; j < m_hessian.size2(); ++j) + { + g_log.information() << std::setw(10); + g_log.information() << m_hessian.get(i,j) << " "; + } + g_log.information() << "\n"; + } + g_log.information().flags(prevState); + } + covar = m_hessian; + covar.invert(); + if(g_log.is(Kernel::Logger::Priority::PRIO_INFORMATION)) + { + g_log.information() << "== Covariance matrix (H^-1) ==\n"; + std::ios::fmtflags prevState = g_log.information().flags(); + g_log.information() << std::left << std::fixed; + for(size_t i = 0; i < covar.size1(); ++i) + { + for(size_t j = 0; j < covar.size2(); ++j) + { + g_log.information() << std::setw(10); + g_log.information() << covar.get(i,j) << " "; + } + g_log.information() << "\n"; + } + g_log.information().flags(prevState); + } + +} + +} // namespace CurveFitting +} // namespace Mantid diff --git a/Code/Mantid/Framework/CurveFitting/src/SeqDomain.cpp b/Code/Mantid/Framework/CurveFitting/src/SeqDomain.cpp index 4a7ab41d8746..b3820a8e73c7 100644 --- a/Code/Mantid/Framework/CurveFitting/src/SeqDomain.cpp +++ b/Code/Mantid/Framework/CurveFitting/src/SeqDomain.cpp @@ -96,7 +96,26 @@ void SeqDomain::leastSquaresVal(const CostFuncLeastSquares& leastSquares) leastSquares.addVal( domain, values ); } } - +/** + * Calculate the value of a simplex cost function + * @param leastSquares :: The least squares cost func to calculate the value for + */ +void SeqDomain::simplexVal(const CostFuncSimplex& simplex) +{ + API::FunctionDomain_sptr domain; + API::FunctionValues_sptr values; + const size_t n = getNDomains(); + for(size_t i = 0; i < n; ++i) + { + values.reset(); + getDomainAndValues( i, domain, values ); + if (!values) + { + throw std::runtime_error("Simplex: undefined FunctionValues."); + } + simplex.addVal( domain, values ); + } +} //------------------------------------------------------------------------------------------------ /** * Calculate the value of a least squares cost function @@ -143,6 +162,30 @@ void SeqDomain::leastSquaresValDerivHessian(const CostFuncLeastSquares& leastSqu } } +/** + * Calculate the value, first and second derivatives of a simplex cost function + * @param simplex :: The simplex cost func to calculate the value for + * @param evalFunction :: Flag to evaluate the value of the cost function + * @param evalDeriv :: Flag to evaluate the first derivatives + * @param evalHessian :: Flag to evaluate the Hessian (second derivatives) + */ +void SeqDomain::simplexValDerivHessian(const CostFuncSimplex& simplex, bool evalFunction, bool evalDeriv, bool evalHessian) +{ + API::FunctionDomain_sptr domain; + API::FunctionValues_sptr values; + const size_t n = getNDomains(); + for(size_t i = 0; i < n; ++i) + { + values.reset(); + getDomainAndValues( i, domain, values ); + if (!values) + { + throw std::runtime_error("Simplex: undefined FunctionValues."); + } + simplex.addValDerivHessian(simplex.getFittingFunction(),domain,values,evalFunction,evalDeriv,evalHessian); + } +} + /** * Calculate the value, first and second derivatives of a RWP cost function * @param rwp :: The rwp cost func to calculate the value for diff --git a/Code/Mantid/Framework/Geometry/inc/MantidGeometry/Crystal/IndexingUtils.h b/Code/Mantid/Framework/Geometry/inc/MantidGeometry/Crystal/IndexingUtils.h index c253550879ad..56428574bf60 100644 --- a/Code/Mantid/Framework/Geometry/inc/MantidGeometry/Crystal/IndexingUtils.h +++ b/Code/Mantid/Framework/Geometry/inc/MantidGeometry/Crystal/IndexingUtils.h @@ -92,7 +92,7 @@ class MANTID_GEOMETRY_DLL IndexingUtils /// Find the UB matrix that most nearly maps hkl to qxyz for 3 or more peaks static double Optimize_UB( Kernel::DblMatrix & UB, const std::vector & hkl_vectors, - const std::vector & q_vectors ); + const std::vector & q_vectors, double* out = NULL); /// Find the vector that best corresponds to plane normal, given 1-D indices static double Optimize_Direction( Kernel::V3D & best_vec, @@ -273,8 +273,6 @@ class MANTID_GEOMETRY_DLL IndexingUtils /// Get a formatted string listing the lattice parameters and cell volume static std::string GetLatticeParameterString( const Kernel::DblMatrix & UB ); - - private: /// Static reference to the logger class diff --git a/Code/Mantid/Framework/Geometry/src/Crystal/IndexingUtils.cpp b/Code/Mantid/Framework/Geometry/src/Crystal/IndexingUtils.cpp index d200788ccabb..017425af1831 100644 --- a/Code/Mantid/Framework/Geometry/src/Crystal/IndexingUtils.cpp +++ b/Code/Mantid/Framework/Geometry/src/Crystal/IndexingUtils.cpp @@ -721,6 +721,7 @@ double IndexingUtils::Optimize_UB( DblMatrix & UB, @param q_vectors std::vector of V3D objects that contains the list of q_vectors that are indexed by the corresponding hkl vectors. + @param out optional array of doubles of residuals for fitting NOTE: The number of hkl_vectors and q_vectors must be the same, and must be at least 3. @@ -739,7 +740,7 @@ double IndexingUtils::Optimize_UB( DblMatrix & UB, */ double IndexingUtils::Optimize_UB( DblMatrix & UB, const std::vector & hkl_vectors, - const std::vector & q_vectors ) + const std::vector & q_vectors , double* out) { if ( UB.numRows() != 3 || UB.numCols() != 3 ) { @@ -815,6 +816,7 @@ double IndexingUtils::Optimize_UB( DblMatrix & UB, for ( size_t i = 0; i < q_vectors.size(); i++ ) { sum_sq_error += gsl_vector_get(residual, i) * gsl_vector_get(residual, i); + if (out != NULL)out[row+i*3] = gsl_vector_get(residual, i); } }