Skip to content

Commit

Permalink
Refs #10199 use same cost function as ISAW
Browse files Browse the repository at this point in the history
  • Loading branch information
VickieLynch committed Nov 11, 2014
1 parent b5d5816 commit 0cd3ae9
Show file tree
Hide file tree
Showing 9 changed files with 746 additions and 25 deletions.
12 changes: 12 additions & 0 deletions Code/Mantid/Framework/Crystal/src/OptimizeCrystalPlacement.cpp
Expand Up @@ -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 );
Expand All @@ -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);
Expand Down
47 changes: 27 additions & 20 deletions Code/Mantid/Framework/Crystal/src/PeakHKLErrors.cpp
Expand Up @@ -21,12 +21,14 @@
#include "MantidKernel/Matrix.h"
#include "MantidKernel/Quat.h"
#include "MantidKernel/V3D.h"
#include "MantidGeometry/Crystal/IndexingUtils.h"
#include <cctype>
#include <string>

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;
Expand Down Expand Up @@ -364,11 +366,7 @@ namespace Mantid

std::map<int, Mantid::Kernel::Matrix<double> > 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");
Expand All @@ -379,6 +377,7 @@ namespace Mantid


double ChiSqTot = 0.0;
std::vector<V3D> q_vectors;
for( size_t i = 0; i< nData; i+= 3 )
{
int peakNum = (int)( .5 + xValues[i] );
Expand All @@ -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<double> sigabc(1);// only need chisq
std::vector<V3D> miller_ind;
std::vector<V3D> 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<double>(num_indexed)*static_cast<double>(nData);
out[p] -= static_cast<double>(num_indexed)/static_cast<double>(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-----------------------------------------------"<<std::endl;
for( size_t p = 0 ; p < nParams() ; p++ )
{
g_log.debug() << parameterName(p)<<"("<<getParameter(p)<<"),";
g_log.debug() << parameterName(p)<<"("<<getParameter(p
)<<"),";
if ((p + 1) % 6 ==0 )
g_log.debug() << std::endl;
}
Expand Down
2 changes: 2 additions & 0 deletions Code/Mantid/Framework/CurveFitting/CMakeLists.txt
Expand Up @@ -22,6 +22,7 @@ set ( SRC_FILES
src/CostFuncFitting.cpp
src/CostFuncLeastSquares.cpp
src/CostFuncRwp.cpp
src/CostFuncSimplex.cpp
src/CubicSpline.cpp
src/DampingMinimizer.cpp
src/DeltaFunction.cpp
Expand Down Expand Up @@ -127,6 +128,7 @@ set ( INC_FILES
inc/MantidCurveFitting/CostFuncFitting.h
inc/MantidCurveFitting/CostFuncLeastSquares.h
inc/MantidCurveFitting/CostFuncRwp.h
inc/MantidCurveFitting/CostFuncSimplex.h
inc/MantidCurveFitting/CubicSpline.h
inc/MantidCurveFitting/DampingMinimizer.h
inc/MantidCurveFitting/DeltaFunction.h
Expand Down
@@ -0,0 +1,115 @@
#ifndef MANTID_CURVEFITTING_COSTFUNCSIMPLEX_H_
#define MANTID_CURVEFITTING_COSTFUNCSIMPLEX_H_

//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidCurveFitting/CostFuncFitting.h"
#include "MantidCurveFitting/GSLMatrix.h"
#include "MantidCurveFitting/GSLVector.h"

namespace Mantid
{
namespace CurveFitting
{
class SeqDomain;
class ParDomain;

/** Cost function for simplex
@author Anders Markvardsen, ISIS, RAL
@date 11/05/2010
Copyright &copy; 2010 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory
This file is part of Mantid.
Mantid is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 3 of the License, or
(at your option) any later version.
Mantid is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
File change history is stored at: <https://github.com/mantidproject/mantid>.
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
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<double>& 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<double>& 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<double> 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_*/
Expand Up @@ -10,6 +10,7 @@
#include "MantidAPI/IDomainCreator.h"
#include "MantidCurveFitting/CostFuncLeastSquares.h"
#include "MantidCurveFitting/CostFuncRwp.h"
#include "MantidCurveFitting/CostFuncSimplex.h"

#include <stdexcept>
#include <vector>
Expand Down Expand Up @@ -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
Expand Down

0 comments on commit 0cd3ae9

Please sign in to comment.