diff --git a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/FABADAMinimizer.h b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/FABADAMinimizer.h
index 96c17435ca3dcc78da315df6bd3ab5dfe0435703..9e389f624b295d9908edbef0822dc4accbf6c9be 100644
--- a/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/FABADAMinimizer.h
+++ b/Code/Mantid/Framework/CurveFitting/inc/MantidCurveFitting/FABADAMinimizer.h
@@ -40,53 +40,52 @@ public:
/// Constructor
FABADAMinimizer();
virtual ~FABADAMinimizer();
- /// Lo ha puesto Roman y no se para que!! creo que es el destructor
-
/// Name of the minimizer.
std::string name() const { return "FABADA"; }
/// Initialize minimizer, i.e. pass a function to minimize.
- virtual void initialize(API::ICostFunction_sptr function,
- size_t maxIterations);
+ virtual void initialize(API::ICostFunction_sptr function,size_t maxIterations);
/// Do one iteration.
virtual bool iterate(size_t iter);
/// Return current value of the cost function
virtual double costFunctionVal();
-private:
- /// Pointer to the cost function. Must be the least squares.
- /// Intentar encontrar una manera de sacar aqui el numero de parametros que
- /// no sea necesaria la cost function
- boost::shared_ptr<CostFuncLeastSquares> m_leastSquares;
- /// The number of iterations done.
- size_t m_counter;
- ///
- size_t m_numberIterations;
- /// The number of changes done in each parameter.
- std::vector<double> m_changes;
- /// The jump for each parameter
- std::vector<double> m_jump;
- /// Parameters.
- GSLVector m_parameters;
- /// Markov chain.
- std::vector<std::vector<double>> m_chain;
- /// The chi square result of previous iteration;
- double m_chi2;
- /// Boolean that indicates if converged
- bool m_converged;
- /// The point when convergence starts
- size_t m_conv_point;
- /// Convergence of each parameter
- std::vector<bool> m_par_converged;
- /// Lower bound for each parameter
- std::vector<double> m_lower;
- /// Upper bound for each parameter
- std::vector<double> m_upper;
- /// Bool that indicates if there is any boundary constraint
- std::vector<bool> m_bound;
- /// Convergence criteria for each parameter
- std::vector<double> m_criteria;
-};
+
+ private:
+ /// Pointer to the cost function. Must be the least squares.
+ /// Intentar encontrar una manera de sacar aqui el numero de parametros que no sea necesaria la cost function
+ boost::shared_ptr<CostFuncLeastSquares> m_leastSquares;
+ /// The number of iterations done.
+ size_t m_counter;
+ ///
+ size_t m_numberIterations;
+ /// The number of changes done in each parameter.
+ std::vector<double> m_changes;
+ /// The jump for each parameter
+ std::vector<double> m_jump;
+ /// Parameters.
+ GSLVector m_parameters;
+ /// Markov chain.
+ std::vector<std::vector<double>> m_chain;
+ /// The chi square result of previous iteration;
+ double m_chi2;
+ /// Boolean that indicates if converged
+ bool m_converged;
+ /// The point when convergence starts
+ size_t m_conv_point;
+ /// Convergence of each parameter
+ std::vector<bool> m_par_converged;
+ ///Lower bound for each parameter
+ std::vector<double> m_lower;
+ ///Upper bound for each parameter
+ std::vector<double> m_upper;
+ /// Bool that indicates if there is any boundary constraint
+ std::vector<bool> m_bound;
+ /// Convergence criteria for each parameter
+ std::vector<double> m_criteria;
+ /// Maximum number of iterations
+ size_t m_max_iter;
+ };
} // namespace CurveFitting
} // namespace Mantid
diff --git a/Code/Mantid/Framework/CurveFitting/src/FABADAMinimizer.cpp b/Code/Mantid/Framework/CurveFitting/src/FABADAMinimizer.cpp
index 96696fb9f968ab9b92dfeb4ce5bcc19e2f97dc4d..1bce0c62f840b7187bb5fb7ab3bda289979f92e7 100644
--- a/Code/Mantid/Framework/CurveFitting/src/FABADAMinimizer.cpp
+++ b/Code/Mantid/Framework/CurveFitting/src/FABADAMinimizer.cpp
@@ -37,8 +37,6 @@ namespace {
Kernel::Logger g_log("FABADAMinimizer");
// absolute maximum number of iterations when fit must converge
const size_t convergenceMaxIterations = 50000;
-// histogram length for the PDF output workspace
-const size_t pdf_length = 50;
// number of iterations when convergence isn't expected
const size_t lowerIterationLimit = 350;
}
@@ -47,37 +45,39 @@ DECLARE_FUNCMINIMIZER(FABADAMinimizer, FABADA)
//----------------------------------------------------------------------------------------------
/// Constructor
-FABADAMinimizer::FABADAMinimizer() : API::IFuncMinimizer(), m_conv_point(0) {
- declareProperty("ChainLength", static_cast<size_t>(10000),
- "Length of the converged chain.");
- declareProperty(
- "ConvergenceCriteria", 0.0001,
- "Variance in Chi square for considering convergence reached.");
- declareProperty(new API::WorkspaceProperty<>("OutputWorkspacePDF", "pdf",
- Kernel::Direction::Output),
- "The name to give the output workspace");
- declareProperty(new API::WorkspaceProperty<>("OutputWorkspaceChain", "chain",
- Kernel::Direction::Output),
- "The name to give the output workspace");
- declareProperty(new API::WorkspaceProperty<>("OutputWorkspaceConverged", "",
- Kernel::Direction::Output,
- API::PropertyMode::Optional),
- "The name to give the output workspace");
- declareProperty(new API::WorkspaceProperty<API::ITableWorkspace>(
- "ChiSquareTable", "chi2", Kernel::Direction::Output),
- "The name to give the output workspace");
- declareProperty(new API::WorkspaceProperty<API::ITableWorkspace>(
- "PdfError", "pdfE", Kernel::Direction::Output),
- "The name to give the output workspace");
-}
+FABADAMinimizer::FABADAMinimizer()
+ {
+ declareProperty("ChainLength",static_cast<size_t>(10000),"Length of the converged chain.");
+ declareProperty("StepsBetweenValues",static_cast<size_t>(10),"Steps realized between keeping each result.");
+ declareProperty("ConvergenceCriteria",0.01,"Variance in Cost Function for considering convergence reached.");
+ declareProperty("JumpAcceptanceRate",0.6666666,"Desired jumping acceptance rate");
+ declareProperty(
+ new API::WorkspaceProperty<>("PDF","PDF",Kernel::Direction::Output),
+ "The name to give the output workspace");
+ declareProperty(
+ new API::WorkspaceProperty<>("Chains","Chain",Kernel::Direction::Output),
+ "The name to give the output workspace");
+ declareProperty(
+ new API::WorkspaceProperty<>("ConvergedChain","ConvergedChain",Kernel::Direction::Output,API::PropertyMode::Optional),
+ "The name to give the output workspace");
+ declareProperty(
+ new API::WorkspaceProperty<API::ITableWorkspace>("CostFunctionTable","CostFunction",Kernel::Direction::Output),
+ "The name to give the output workspace");
+ declareProperty(
+ new API::WorkspaceProperty<API::ITableWorkspace>("Parameters","Parameters",Kernel::Direction::Output),
+ "The name to give the output workspace");
+
+
+ }
+
+
+ //----------------------------------------------------------------------------------------------
-//----------------------------------------------------------------------------------------------
/// Destructor
FABADAMinimizer::~FABADAMinimizer() {}
/// Initialize minimizer. Set initial values for all private members.
-void FABADAMinimizer::initialize(API::ICostFunction_sptr function,
- size_t maxIterations) {
+void FABADAMinimizer::initialize(API::ICostFunction_sptr function,size_t maxIterations) {
m_leastSquares = boost::dynamic_pointer_cast<CostFuncLeastSquares>(function);
if (!m_leastSquares) {
@@ -89,10 +89,12 @@ void FABADAMinimizer::initialize(API::ICostFunction_sptr function,
m_leastSquares->getParameters(m_parameters);
API::IFunction_sptr fun = m_leastSquares->getFittingFunction();
+
if (fun->nParams() == 0) {
throw std::invalid_argument("Function has 0 fitting parameters.");
}
+
size_t n = getProperty("ChainLength");
m_numberIterations = n / fun->nParams();
@@ -131,9 +133,14 @@ void FABADAMinimizer::initialize(API::ICostFunction_sptr function,
}
}
}
+ else {
+ m_lower.push_back(-10e100);
+ m_upper.push_back(10e100);
+ }
std::vector<double> v;
v.push_back(p);
m_chain.push_back(v);
+ m_max_iter = maxIterations;
m_changes.push_back(0);
m_par_converged.push_back(false);
m_criteria.push_back(getProperty("ConvergenceCriteria"));
@@ -148,6 +155,7 @@ void FABADAMinimizer::initialize(API::ICostFunction_sptr function,
v.push_back(m_chi2);
m_chain.push_back(v);
m_converged = false;
+ m_max_iter = maxIterations;
}
/// Do one iteration. Returns true if iterations to be continued, false if they
@@ -161,52 +169,59 @@ bool FABADAMinimizer::iterate(size_t) {
size_t n = m_leastSquares->nParams();
size_t m = n;
- // Just for the last iteration. For doing exactly the indicated number of
- // iterations.
- if (m_converged && m_counter == (m_numberIterations)) {
- size_t t = getProperty("ChainLength");
- m = t % n;
- }
-
- // Do one iteration of FABADA's algorithm for each parameter.
- for (size_t i = 0; i < m; i++) {
- GSLVector new_parameters = m_parameters;
-
- // Calculate the step, depending on convergence reached or not
- double step;
- if (m_converged) {
- boost::mt19937 mt;
- mt.seed(123 * (int(m_counter) + 45 * int(i)));
-
- boost::normal_distribution<double> distr(0.0, std::abs(m_jump[i]));
- boost::variate_generator<
- boost::mt19937, boost::normal_distribution<double>> gen(mt, distr);
-
- step = gen();
-
- } else {
- step = m_jump[i];
+ // Just for the last iteration. For doing exactly the indicated number of iterations.
+ if(m_converged && m_counter == m_numberIterations){
+ size_t t = getProperty("ChainLength");
+ m = t % n;
}
- // Couts just for helping when developing when coding
- /// std::cout << std::endl << m_counter << "." << i <<
- /// std::endl<<std::endl<< m_parameters.get(i)<<std::endl; //DELETE AT THE
- /// END
- /// std::cout << "Real step: " << step << " Due to the m_jump: " <<
- /// m_jump[i] << std::endl; //DELETE AT THE END
-
- // Calculate the new value of the parameter
- double new_value = m_parameters.get(i) + step;
-
- // Comproves if it is inside the boundary constrinctions. If not, changes
- // it.
- if (m_bound[i]) {
- if (new_value < m_lower[i]) {
- new_value = m_lower[i] + (m_lower[i] - new_value) / 2;
+
+ // Do one iteration of FABADA's algorithm for each parameter.
+ for(size_t i = 0; i<m ; i++)
+ {
+ GSLVector new_parameters = m_parameters;
+
+ // Calculate the step, depending on convergence reached or not
+ double step;
+ if (m_converged || m_bound[i])
+ {
+ boost::mt19937 mt;
+ mt.seed(123*(int(m_counter)+45*int(i))); //Numeros inventados para la seed
+ boost::random::normal_distribution<> distr(0,std::abs(m_jump[i]));
+ step = distr(mt);
+
}
- if (new_value > m_upper[i]) {
- new_value = m_upper[i] - (new_value - m_upper[i]) / 2;
+ else {
+ step = m_jump[i];
+ }
+
+ // Calculate the new value of the parameter
+ double new_value = m_parameters.get(i) + step;
+
+ // Comproves if it is inside the boundary constrinctions. If not, changes it.
+ if (m_bound[i])
+ {
+ while (new_value < m_lower[i]) {
+ if (std::abs(step) > m_upper[i]-m_lower[i]) {
+ new_value = m_parameters.get(i) + step/10.0;
+ step=step/10;
+ m_jump[i] = m_jump[i]/10;
+ }
+ else {
+ new_value = m_lower[i] + std::abs(step) - (m_parameters.get(i)-m_lower[i]);
+ }
+ }
+ while (new_value > m_upper[i]) {
+ if (std::abs(step) > m_upper[i]-m_lower[i]) {
+ new_value = m_parameters.get(i) + step/10.0;
+ step=step/10;
+ m_jump[i] = m_jump[i]/10;
+ }
+ else {
+ new_value = m_upper[i] - (std::abs(step) + m_parameters.get(i)-m_upper[i]);
+ }
+ }
}
- }
+
// Set the new value in order to calculate the new Chi square value
if (boost::math::isnan(new_value)) {
@@ -216,8 +231,6 @@ bool FABADAMinimizer::iterate(size_t) {
m_leastSquares->setParameter(i, new_value);
double chi2_new = m_leastSquares->val();
- /// std::cout << "OLD Chi2: " << m_chi2 << " NEW Chi2: " << chi2_new
- /// << std::endl; // DELETE AT THE END
// If new Chi square value is lower, jumping directly to new parameter
if (chi2_new < m_chi2) {
@@ -228,22 +241,19 @@ bool FABADAMinimizer::iterate(size_t) {
m_parameters = new_parameters;
m_chi2 = chi2_new;
m_changes[i] += 1;
- /// std::cout << "Salta directamente!!" << std::endl;// DELETE AT THE END
+
}
// If new Chi square value is higher, it depends on the probability
else {
// Calculate probability of change
double prob = exp((m_chi2 / 2.0) - (chi2_new / 2.0));
- /// std::cout << "PROBABILIDAD cambio: " << prob << std::endl;// DELETE
- /// AT THE END
// Decide if changing or not
boost::mt19937 mt;
mt.seed(int(time_t()) + 48 * (int(m_counter) + 76 * int(i)));
boost::uniform_real<> distr(0.0, 1.0);
double p = distr(mt);
- /// std::cout << " Random number " << p << std::endl;// DELETE AT THE END
if (p <= prob) {
for (size_t j = 0; j < n; j++) {
m_chain[j].push_back(new_parameters.get(j));
@@ -252,7 +262,6 @@ bool FABADAMinimizer::iterate(size_t) {
m_parameters = new_parameters;
m_chi2 = chi2_new;
m_changes[i] += 1;
- /// std::cout << "SI hay cambio" << std::endl;// DELETE AT THE END
} else {
for (size_t j = 0; j < n; j++) {
m_chain[j].push_back(m_parameters.get(j));
@@ -260,91 +269,71 @@ bool FABADAMinimizer::iterate(size_t) {
m_chain[n].push_back(m_chi2);
m_leastSquares->setParameter(i, new_value - m_jump[i]);
m_jump[i] = -m_jump[i];
- /// std::cout << "NO hay cambio" << std::endl;// DELETE AT THE END
}
}
- /// std::cout << std::endl << std::endl << std::endl;// DELETE AT THE END
-
- const size_t jumpCheckingRate = 200;
- const double lowJumpLimit = 1e-15;
-
- // Update the jump once each jumpCheckingRate iterations
- if (m_counter % jumpCheckingRate == 150) {
- double jnew;
- // All this is just a temporal test...
- std::vector<double>::const_iterator first = m_chain[n].end() - 41;
- std::vector<double>::const_iterator last = m_chain[n].end();
- std::vector<double> test(first, last);
- int c = 0;
- for (int j = 0; j < 39; ++j) {
- if (test[j] == test[j + 1]) {
- c += 1;
- }
- }
- if (c > 38) {
- jnew = m_jump[i] / 100;
- } // ...untill here.
- else {
- if (m_changes[i] == 0.0) {
- jnew = m_jump[i] / 10.0;
- } else {
- double f = m_changes[i] / double(m_counter);
- jnew = m_jump[i] * f / 0.6666666666;
- /// std::cout << f << " m_counter "<< m_counter << " m_changes
- /// " << m_changes[i] << std::endl; // DELETE AT THE END
- }
- }
- m_jump[i] = jnew;
+ const size_t jumpCheckingRate = 200;
+ const double lowJumpLimit = 1e-25;
+ const double jumpAR = getProperty("JumpAcceptanceRate");
- // Check if the new jump is too small. It means that it has been a wrong
- // convergence.
- if (std::abs(m_jump[i]) < lowJumpLimit) {
- API::IFunction_sptr fun = m_leastSquares->getFittingFunction();
- throw std::runtime_error(
- "Wrong convergence for parameter " + fun->parameterName(i) +
- ". Try to set a proper initial value this parameter");
+ // Update the jump once each jumpCheckingRate iterations
+ if (m_counter % jumpCheckingRate == 150) //JUMP CHECKING RATE IS 200, BUT IS NOT CHECKED AT FIRST STEP, IT IS AT 150
+ {
+ double jnew;
+ if (m_changes[i] == 0.0) {
+ jnew = m_jump[i]/10.0; //JUST FOR THE CASE THERE HAS NOT BEEN ANY CHANGE.
+ }
+ else {
+ double f = m_changes[i]/double(m_counter);
+ jnew = m_jump[i]*f/jumpAR;
+ }
+
+ m_jump[i] = jnew;
+
+ // Check if the new jump is too small. It means that it has been a wrong convergence.
+ if (std::abs(m_jump[i])<lowJumpLimit) {
+ API::IFunction_sptr fun = m_leastSquares->getFittingFunction();
+ g_log.warning()<< "Wrong convergence for parameter " + fun -> parameterName(i) +". Try to set a proper initial value for this parameter"<< std::endl;
+ }
}
- }
- // Check if the Chi square value has converged for parameter i.
- if (!m_par_converged[i] &&
- m_counter > 350) // It only check since the iteration number 350
- {
- if (chi2_new != m_chi2) {
- double chi2_quotient = std::abs(chi2_new - m_chi2) / m_chi2;
- if (chi2_quotient < m_criteria[i]) {
- m_par_converged[i] = true;
+ // Check if the Chi square value has converged for parameter i.
+ const size_t startingPoint = 350; // The iteration since it starts to check if convergence is reached
+ if (!m_par_converged[i] && m_counter > startingPoint)
+ {
+ if (chi2_new != m_chi2) {
+ double chi2_quotient = std::abs(chi2_new-m_chi2)/m_chi2;
+ if (chi2_quotient < m_criteria[i]){
+ m_par_converged[i] = true;
+ }
}
}
}
- }
- m_counter +=
- 1; // Update the counter, after finishing the iteration for each parameter
+ m_counter += 1; // Update the counter, after finishing the iteration for each parameter
- // Check if Chi square has converged for all the parameters.
- if (m_counter > lowerIterationLimit && !m_converged) {
- size_t t = 0;
- for (size_t i = 0; i < n; i++) {
- if (m_par_converged[i]) {
- t += 1;
+ // Check if Chi square has converged for all the parameters.
+ if (m_counter > lowerIterationLimit && !m_converged) {
+ size_t t = 0;
+ for (size_t i = 0; i < n; i++) {
+ if (m_par_converged[i]) {
+ t += 1;
+ }
}
- }
- // If all parameters have converged:
- // It set up both the counter and the changes' vector to 0, in order to
- // consider only the
- // data of the converged part of the chain, when updating the jump.
- if (t == n) {
- m_converged = true;
- m_conv_point = m_counter * n + 1;
- m_counter = 0;
- for (size_t i = 0; i < n; ++i) {
- m_changes[i] = 0;
+ // If all parameters have converged:
+ // It set up both the counter and the changes' vector to 0, in order to
+ // consider only the data of the converged part of the chain, when updating the jump.
+ if (t == n) {
+ m_converged = true;
+ m_conv_point = m_counter * n + 1;
+ m_counter = 0;
+ for (size_t i = 0; i < n; ++i) {
+ m_changes[i] = 0;
+ }
}
}
- }
+
if (!m_converged) {
// If there is not convergence continue the iterations.
@@ -365,186 +354,266 @@ bool FABADAMinimizer::iterate(size_t) {
failed.replace(failed.end() - 2, failed.end(), ".");
throw std::runtime_error(
"Convegence NOT reached after " +
- boost::lexical_cast<std::string>(m_counter) +
+ boost::lexical_cast<std::string>(m_max_iter) +
" iterations.\n Try to set better initial values for parameters: " +
failed);
}
- } else {
+ }
+
+ else {
// If convergence has been reached, continue untill complete the chain
// length.
if (m_counter <= m_numberIterations) {
return true;
}
- // When the all the iterations have been done, calculate and show all the
- // results.
- else {
- // Create the workspace for the Probability Density Functions
- API::MatrixWorkspace_sptr ws = API::WorkspaceFactory::Instance().create(
- "Workspace2D", n, pdf_length + 1, pdf_length);
-
- // Create the workspace for the parameters' value and errors.
- API::ITableWorkspace_sptr wsPdfE =
- API::WorkspaceFactory::Instance().createTable("TableWorkspace");
- wsPdfE->addColumn("str", "Name");
- wsPdfE->addColumn("double", "Value");
- wsPdfE->addColumn("double", "Left's error");
- wsPdfE->addColumn("double", "Rigth's error");
-
- std::vector<double>::iterator pos_min = std::min_element(
- m_chain[n].begin() + m_conv_point,
- m_chain[n]
- .end()); // Calculate the position of the minimum Chi aquare value
- m_chi2 = *pos_min;
- // index of the minimum chi^2
- size_t minIndex =
- static_cast<size_t>(std::distance(m_chain[n].begin(), pos_min));
-
- std::vector<double> par_def(n);
- API::IFunction_sptr fun = m_leastSquares->getFittingFunction();
+ //If convergence has been reached, but the maximum of iterations have been reached before finishing the chain, stop and throw the error.
+ if (m_counter >= m_max_iter)
+ {
+ throw std::length_error("Convegence reached but Max Iterations parameter insufficient for creating the whole chain.\n Increase Max Iterations");
+ return false;
+ }
+
+ // When the all the iterations have been done, calculate and show all the results.
+ else
+ {
- // Do one iteration for each parameter.
- for (size_t j = 0; j < n; ++j) {
- // Calculate the parameter value and the errors
- par_def[j] = m_chain[j][minIndex];
- std::vector<double>::const_iterator first =
- m_chain[j].begin() + m_conv_point;
- std::vector<double>::const_iterator last = m_chain[j].end();
+ //Creating the reduced chain (considering only one each "Steps between values" values)
+ size_t cl = getProperty("ChainLength");
+ size_t n_steps = getProperty("StepsBetweenValues");
+ size_t conv_length = size_t(double(cl)/double(n_steps));
+ std::vector<std::vector<double>> red_conv_chain;
+ for(size_t e=0;e<=n;++e) {
+ std::vector<double> v;
+ v.push_back(m_chain[e][m_conv_point]);
+ red_conv_chain.push_back(v);
+ }
+
+ //Calculate the red_conv_chain for the cost fuction.
+ std::vector<double>::const_iterator first = m_chain[n].begin()+m_conv_point;
+ std::vector<double>::const_iterator last = m_chain[n].end();
std::vector<double> conv_chain(first, last);
- size_t conv_length = conv_chain.size();
- std::sort(conv_chain.begin(), conv_chain.end());
- std::vector<double>::const_iterator pos_par =
- std::find(conv_chain.begin(), conv_chain.end(), par_def[j]);
- int sigma = int(0.34 * double(conv_length));
- // make sure the iterator is valid in any case
- std::vector<double>::const_iterator pos_left = conv_chain.begin();
- if (sigma < static_cast<int>(std::distance(pos_left, pos_par))) {
- pos_left = pos_par - sigma;
+ for(size_t k=1; k<conv_length; ++k) {
+ red_conv_chain[n].push_back(conv_chain[n_steps*k]);
}
- // make sure the iterator is valid in any case
- std::vector<double>::const_iterator pos_right = conv_chain.end() - 1;
- if (sigma < static_cast<int>(std::distance(pos_par, pos_right))) {
- pos_right = pos_par + sigma;
+
+ std::vector<double>::iterator pos_min = std::min_element(red_conv_chain[n].begin(),red_conv_chain[n].end()); // Calculate the position of the minimum Chi square value
+ m_chi2 = *pos_min;
+
+
+ std::vector<double> par_def(n);
+ std::vector<double> error_left(n);
+ std::vector<double> error_rigth(n);
+ API::IFunction_sptr fun = m_leastSquares->getFittingFunction();
+
+
+ // Do one iteration for each parameter.
+ for (size_t j =0; j < n; ++j)
+ {
+ // Calculate the parameter value and the errors
+ std::vector<double>::const_iterator first = m_chain[j].begin()+m_conv_point;
+ std::vector<double>::const_iterator last = m_chain[j].end();
+ std::vector<double> conv_chain(first, last);
+ for(size_t k=0; k<conv_length; ++k) {
+ red_conv_chain[j].push_back(conv_chain[n_steps*k]);
+ }
+ par_def[j]=red_conv_chain[j][pos_min-red_conv_chain[n].begin()];
+ std::sort(red_conv_chain[j].begin(),red_conv_chain[j].end());
+ std::vector<double>::const_iterator pos_par = std::find(red_conv_chain[j].begin(),red_conv_chain[j].end(),par_def[j]);
+ int sigma = int(0.34*double(conv_length));
+ std::vector<double>::const_iterator pos_left = pos_par - sigma;
+ std::vector<double>::const_iterator pos_right = pos_par + sigma;
+ error_left[j]= *pos_left - *pos_par;
+ error_rigth[j]= *pos_right - *pos_par;
+
}
- API::TableRow row = wsPdfE->appendRow();
- row << fun->parameterName(j) << par_def[j] << *pos_left - *pos_par
- << *pos_right - *pos_par;
-
- // Calculate the Probability Density Function
- std::vector<double> pdf_y(pdf_length, 0);
- double start = conv_chain[0];
- double bin = (conv_chain[conv_length - 1] - start) / pdf_length;
- size_t step = 0;
- MantidVec &X = ws->dataX(j);
- MantidVec &Y = ws->dataY(j);
- X[0] = start;
- for (size_t i = 1; i < pdf_length + 1; i++) {
- double bin_end = start + static_cast<double>(i) * bin;
- X[i] = bin_end;
- while (step < conv_length && conv_chain[step] <= bin_end) {
- pdf_y[i - 1] += 1;
- ++step;
- }
- Y[i - 1] = pdf_y[i - 1] / (double(conv_length) * bin);
+
+ const bool cond1 = !getPropertyValue("Parameters").empty();
+
+ if (cond1) {
+
+ // Create the workspace for the parameters' value and errors.
+ API::ITableWorkspace_sptr wsPdfE = API::WorkspaceFactory::Instance().createTable("TableWorkspace");
+ wsPdfE -> addColumn("str","Name");
+ wsPdfE -> addColumn("double","Value");
+ wsPdfE -> addColumn("double","Left's error");
+ wsPdfE -> addColumn("double","Rigth's error");
+
+ for (size_t j =0; j < n; ++j) {
+ API::TableRow row = wsPdfE -> appendRow();
+ row << fun->parameterName(j) << par_def[j] << error_left[j] << error_rigth[j];
+ }
+ // Set and name the Parameter Errors workspace.
+ setProperty("Parameters",wsPdfE);
}
- // Calculate the most probable value, from the PDF.
- std::vector<double>::iterator pos_MP =
- std::max_element(pdf_y.begin(), pdf_y.end());
- double mostP = X[pos_MP - pdf_y.begin()] + (bin / 2.0);
- m_leastSquares->setParameter(j, mostP);
- }
- // Set and name the two workspaces already calculated.
- setProperty("OutputWorkspacePDF", ws);
- setProperty("PdfError", wsPdfE);
-
- // Create the workspace for the complete parameters' chain (the last
- // histogram is for the Chi square).
- size_t chain_length = m_chain[0].size();
- API::MatrixWorkspace_sptr wsC = API::WorkspaceFactory::Instance().create(
- "Workspace2D", n + 1, chain_length, chain_length);
-
- // Do one iteration for each parameter plus one for Chi square.
- for (size_t j = 0; j < n + 1; ++j) {
- MantidVec &X = wsC->dataX(j);
- MantidVec &Y = wsC->dataY(j);
- for (size_t k = 0; k < chain_length; ++k) {
- X[k] = double(k);
- Y[k] = m_chain[j][k];
+ // Set the best parameter values
+ for (size_t j =0; j<n; ++j){
+ m_leastSquares -> setParameter(j,par_def[j]);
}
- }
- // Set and name the workspace for the complete chain
- setProperty("OutputWorkspaceChain", wsC);
-
- // Read if necessary to show the workspace for the converged part of the
- // chain.
- const bool con = !getPropertyValue("OutputWorkspaceConverged").empty();
-
- if (con) {
- // Create the workspace for the converged part of the chain.
- size_t conv_length = (m_counter - 1) * n + m;
- API::MatrixWorkspace_sptr wsConv =
- API::WorkspaceFactory::Instance().create("Workspace2D", n + 1,
- conv_length, conv_length);
-
- // Do one iteration for each parameter plus one for Chi square.
- for (size_t j = 0; j < n + 1; ++j) {
- std::vector<double>::const_iterator first =
- m_chain[j].begin() + m_conv_point;
- std::vector<double>::const_iterator last = m_chain[j].end();
- std::vector<double> conv_chain(first, last);
- MantidVec &X = wsConv->dataX(j);
- MantidVec &Y = wsConv->dataY(j);
- for (size_t k = 0; k < conv_length; ++k) {
- X[k] = double(k);
- Y[k] = conv_chain[k];
+ double mostPchi2;
+ //const bool cond2 = !getPropertyValue("PDF").empty();
+ bool cond2 = true; //This workspace is necessary to be calculated
+
+ if (cond2) {
+
+
+ // Create the workspace for the Probability Density Functions
+ size_t pdf_length = 20; // histogram length for the PDF output workspace
+ API::MatrixWorkspace_sptr ws = API::WorkspaceFactory::Instance().create("Workspace2D",n+1, pdf_length + 1, pdf_length);
+
+ // Calculate the cost function Probability Density Function
+ std::sort(red_conv_chain[n].begin(),red_conv_chain[n].end());
+ std::vector<double> pdf_y(pdf_length,0);
+ double start = red_conv_chain[n][0];
+ double bin = (red_conv_chain[n][conv_length-1] - start)/double(pdf_length);
+ size_t step = 0;
+ MantidVec & X = ws->dataX(n);
+ MantidVec & Y = ws->dataY(n);
+ X[0] = start;
+ for (size_t i = 1; i<pdf_length+1; i++)
+ {
+ double bin_end = start + double(i)*bin;
+ X[i] = bin_end;
+ while(step<conv_length && red_conv_chain[n][step] <= bin_end) {
+ pdf_y[i-1] += 1;
+ ++step;
+ }
+ Y[i-1] = pdf_y[i-1]/(double(conv_length)*bin);
+ }
+
+ std::vector<double>::iterator pos_MPchi2 = std::max_element(pdf_y.begin(),pdf_y.end());
+
+ if(pos_MPchi2-pdf_y.begin() == 0) {
+ //mostPchi2 = X[pos_MPchi2-pdf_y.begin()];
+ mostPchi2 = *pos_min;
+ }
+ else {
+ mostPchi2 = X[pos_MPchi2-pdf_y.begin()]+(bin/2.0);
+ }
+
+ // Do one iteration for each parameter.
+ for (size_t j =0; j < n; ++j) {
+ // Calculate the Probability Density Function
+ std::vector<double> pdf_y(pdf_length,0);
+ double start = red_conv_chain[j][0];
+ double bin = (red_conv_chain[j][conv_length-1] - start)/double(pdf_length);
+ size_t step = 0;
+ MantidVec & X = ws->dataX(j);
+ MantidVec & Y = ws->dataY(j);
+ X[0] = start;
+ for (size_t i = 1; i<pdf_length+1; i++)
+ {
+ double bin_end = start + double(i)*bin;
+ X[i] = bin_end;
+ while(step<conv_length && red_conv_chain[j][step] <= bin_end) {
+ pdf_y[i-1] += 1;
+ ++step;
+ }
+ Y[i-1] = pdf_y[i-1]/(double(conv_length)*bin);
+ }
+
+ // Calculate the most probable value, from the PDF.
+ std::vector<double>::iterator pos_MP = std::max_element(pdf_y.begin(),pdf_y.end());
+ double mostP = X[pos_MP-pdf_y.begin()]+(bin/2.0);
+ m_leastSquares -> setParameter(j,mostP);
+
+ }
+
+ // Set and name the PDF workspace.
+ setProperty("PDF",ws);
+ }
+
+ const bool cond3 = !getPropertyValue("Chains").empty();
+
+ if (cond3) {
+
+ // Create the workspace for the complete parameters' chain (the last histogram is for the Chi square).
+ size_t chain_length = m_chain[0].size();
+ API::MatrixWorkspace_sptr wsC = API::WorkspaceFactory::Instance().create("Workspace2D",n+1, chain_length, chain_length);
+
+ // Do one iteration for each parameter plus one for Chi square.
+ for (size_t j =0; j < n+1; ++j)
+ {
+ MantidVec & X = wsC->dataX(j);
+ MantidVec & Y = wsC->dataY(j);
+ for(size_t k=0; k < chain_length; ++k) {
+ X[k] = double(k);
+ Y[k] = m_chain[j][k];
+ }
+ }
+
+ // Set and name the workspace for the complete chain
+ setProperty("Chains",wsC);
+ }
+
+ // Read if necessary to show the workspace for the converged part of the chain.
+ const bool cond4 = !getPropertyValue("ConvergedChain").empty();
+
+ if (cond4) {
+ // Create the workspace for the converged part of the chain.
+ API::MatrixWorkspace_sptr wsConv = API::WorkspaceFactory::Instance().create("Workspace2D",n+1, conv_length, conv_length);
+
+ // Do one iteration for each parameter plus one for Chi square.
+ for (size_t j =0; j < n+1; ++j)
+ {
+ std::vector<double>::const_iterator first = m_chain[j].begin()+m_conv_point;
+ std::vector<double>::const_iterator last = m_chain[j].end();
+ std::vector<double> conv_chain(first, last);
+ MantidVec & X = wsConv->dataX(j);
+ MantidVec & Y = wsConv->dataY(j);
+ for(size_t k=0; k < conv_length; ++k) {
+ X[k] = double(k);
+ Y[k] = conv_chain[n_steps*k];
+ }
}
+
+ // Set and name the workspace for the converged part of the chain.
+ setProperty("ConvergedChain",wsConv);
}
- // Set and name the workspace for the converged part of the chain.
- setProperty("OutputWorkspaceConverged", wsConv);
- }
+ // Read if necessary to show the workspace for the Chi square values.
+ const bool cond5 = !getPropertyValue("CostFunctionTable").empty();
+
+ if (cond5) {
+
+ // Create the workspace for the Chi square values.
+ API::ITableWorkspace_sptr wsChi2 = API::WorkspaceFactory::Instance().createTable("TableWorkspace");
+ wsChi2 -> addColumn("double","Chi2min");
+ wsChi2 -> addColumn("double","Chi2MP");
+ wsChi2 -> addColumn("double","Chi2min_red");
+ wsChi2 -> addColumn("double","Chi2MP_red");
- // Create the workspace for the Chi square values.
- API::ITableWorkspace_sptr wsChi2 =
- API::WorkspaceFactory::Instance().createTable("TableWorkspace");
- wsChi2->addColumn("double", "Chi2min");
- wsChi2->addColumn("double", "Chi2MP");
- wsChi2->addColumn("double", "Chi2min_red");
- wsChi2->addColumn("double", "Chi2MP_red");
-
- // Calculate de Chi square most probable.
- double Chi2MP = m_leastSquares->val();
-
- // Reset the best parameter values ---> Si al final no se muestra la tabla
- // que sale por defecto, esto se podra borrar...
- for (size_t j = 0; j < n; ++j) {
- m_leastSquares->setParameter(j, par_def[j]);
- }
- // Obtain the quantity of the initial data.
- API::FunctionDomain_sptr domain = m_leastSquares->getDomain();
- size_t data_number = domain->size();
- // Calculate the value for the reduced Chi square.
- double Chi2min_red =
- *pos_min / (double(data_number - n)); // For de minimum value.
- double Chi2MP_red =
- Chi2MP / (double(data_number - n)); // For the most probable.
+ // Obtain the quantity of the initial data.
+ API::FunctionDomain_sptr domain = m_leastSquares -> getDomain();
+ size_t data_number = domain->size();
- // Add the information to the workspace and name it.
- API::TableRow row = wsChi2->appendRow();
- row << *pos_min << Chi2MP << Chi2min_red << Chi2MP_red;
- setProperty("ChiSquareTable", wsChi2);
+ // Calculate the value for the reduced Chi square.
+ double Chi2min_red = m_chi2/(double(data_number-n)); // For de minimum value.
+ double mostPchi2_red = mostPchi2/(double(data_number-n));
- return false;
+ // Add the information to the workspace and name it.
+ API::TableRow row = wsChi2 -> appendRow();
+ row << m_chi2 << mostPchi2 << Chi2min_red << mostPchi2_red;
+ setProperty("CostFunctionTable",wsChi2);
+ }
+
+ // Set the best parameter values
+ for (size_t j =0; j<n; ++j){
+ m_leastSquares -> setParameter(j,par_def[j]);
+ }
+
+ return false;
}
}
return true;
}
-double FABADAMinimizer::costFunctionVal() { return m_chi2; }
+
+ double FABADAMinimizer::costFunctionVal() { return m_chi2; }
} // namespace CurveFitting
} // namespace Mantid
diff --git a/Code/Mantid/Framework/CurveFitting/test/FABADAMinimizerTest.h b/Code/Mantid/Framework/CurveFitting/test/FABADAMinimizerTest.h
index a8f1857f0f065c6b8cbb6b988828c70a02a36323..5e8d39d8735778fa8484c05257288f22b64e3910 100644
--- a/Code/Mantid/Framework/CurveFitting/test/FABADAMinimizerTest.h
+++ b/Code/Mantid/Framework/CurveFitting/test/FABADAMinimizerTest.h
@@ -46,7 +46,7 @@ public:
fit.setProperty("WorkspaceIndex",0);
fit.setProperty("CreateOutput",true);
fit.setProperty("MaxIterations",100000);
- fit.setProperty("Minimizer", "FABADA,ChainLength=5000,ConvergenceCriteria = 0.1, OutputWorkspaceConverged=conv");
+ fit.setProperty("Minimizer", "FABADA,ChainLength=5000,StepsBetweenValues=10,ConvergenceCriteria = 0.1");
TS_ASSERT_THROWS_NOTHING( fit.execute() );
@@ -59,40 +59,40 @@ public:
size_t n = fun -> nParams();
- TS_ASSERT( AnalysisDataService::Instance().doesExist("pdf") );
+ TS_ASSERT( AnalysisDataService::Instance().doesExist("PDF") );
MatrixWorkspace_sptr wsPDF = boost::dynamic_pointer_cast<MatrixWorkspace>(
- API::AnalysisDataService::Instance().retrieve("pdf"));
+ API::AnalysisDataService::Instance().retrieve("PDF"));
TS_ASSERT(wsPDF);
TS_ASSERT_EQUALS(wsPDF->getNumberHistograms(),n);
const Mantid::MantidVec& X = wsPDF->dataX(0);
const Mantid::MantidVec& Y = wsPDF->dataY(0);
- TS_ASSERT_EQUALS(X.size(), 51);
- TS_ASSERT_EQUALS(Y.size(), 50);
-
- TS_ASSERT( AnalysisDataService::Instance().doesExist("chi2") );
- ITableWorkspace_sptr chi2table = boost::dynamic_pointer_cast<ITableWorkspace>(
- API::AnalysisDataService::Instance().retrieve("chi2"));
-
- TS_ASSERT(chi2table);
- TS_ASSERT_EQUALS(chi2table->columnCount(), 4);
- TS_ASSERT_EQUALS(chi2table->rowCount(), 1);
- TS_ASSERT_EQUALS(chi2table->getColumn(0)->type(), "double");
- TS_ASSERT_EQUALS(chi2table->getColumn(0)->name(), "Chi2min");
- TS_ASSERT_EQUALS(chi2table->getColumn(1)->type(), "double");
- TS_ASSERT_EQUALS(chi2table->getColumn(1)->name(), "Chi2MP");
- TS_ASSERT_EQUALS(chi2table->getColumn(2)->type(), "double");
- TS_ASSERT_EQUALS(chi2table->getColumn(2)->name(), "Chi2min_red");
- TS_ASSERT_EQUALS(chi2table->getColumn(3)->type(), "double");
- TS_ASSERT_EQUALS(chi2table->getColumn(3)->name(), "Chi2MP_red");
- TS_ASSERT(chi2table->Double(0,0) <= chi2table->Double(0,1));
- TS_ASSERT(chi2table->Double(0,2) <= chi2table->Double(0,3));
- TS_ASSERT_DELTA(chi2table->Double(0,0), chi2table->Double(0,1), 1);
- TS_ASSERT_DELTA(chi2table->Double(0,0), 0.0, 1.0);
-
- TS_ASSERT( AnalysisDataService::Instance().doesExist("conv") );
+ TS_ASSERT_EQUALS(X.size(), 21);
+ TS_ASSERT_EQUALS(Y.size(), 20);
+
+ TS_ASSERT( AnalysisDataService::Instance().doesExist("CostFunction") );
+ ITableWorkspace_sptr CostFunctionTable = boost::dynamic_pointer_cast<ITableWorkspace>(
+ API::AnalysisDataService::Instance().retrieve("CostFunction"));
+
+ TS_ASSERT(CostFunctionTable);
+ TS_ASSERT_EQUALS(CostFunctionTable->columnCount(), 4);
+ TS_ASSERT_EQUALS(CostFunctionTable->rowCount(), 1);
+ TS_ASSERT_EQUALS(CostFunctionTable->getColumn(0)->type(), "double");
+ TS_ASSERT_EQUALS(CostFunctionTable->getColumn(0)->name(), "Chi2min");
+ TS_ASSERT_EQUALS(CostFunctionTable->getColumn(1)->type(), "double");
+ TS_ASSERT_EQUALS(CostFunctionTable->getColumn(1)->name(), "Chi2MP");
+ TS_ASSERT_EQUALS(CostFunctionTable->getColumn(2)->type(), "double");
+ TS_ASSERT_EQUALS(CostFunctionTable->getColumn(2)->name(), "Chi2min_red");
+ TS_ASSERT_EQUALS(CostFunctionTable->getColumn(3)->type(), "double");
+ TS_ASSERT_EQUALS(CostFunctionTable->getColumn(3)->name(), "Chi2MP_red");
+ TS_ASSERT(CostFunctionTable->Double(0,0) <= CostFunctionTable->Double(0,1));
+ TS_ASSERT(CostFunctionTable->Double(0,2) <= CostFunctionTable->Double(0,3));
+ TS_ASSERT_DELTA(CostFunctionTable->Double(0,0), CostFunctionTable->Double(0,1), 1);
+ TS_ASSERT_DELTA(CostFunctionTable->Double(0,0), 0.0, 1.0);
+
+ TS_ASSERT( AnalysisDataService::Instance().doesExist("ConvergedChain") );
MatrixWorkspace_sptr wsConv = boost::dynamic_pointer_cast<MatrixWorkspace>(
- API::AnalysisDataService::Instance().retrieve("conv"));
+ API::AnalysisDataService::Instance().retrieve("ConvergedChain"));
TS_ASSERT(wsConv);
TS_ASSERT_EQUALS(wsConv->getNumberHistograms(),n+1);
@@ -107,28 +107,27 @@ public:
TS_ASSERT_EQUALS(wsChain->getNumberHistograms(),n+1);
const Mantid::MantidVec& Xchain = wsChain->dataX(0);
- TS_ASSERT_EQUALS(Xchain.size(), 6881);
TS_ASSERT_EQUALS(Xchain[5000], 5000);
TS_ASSERT(Xconv.size() < Xchain.size());
- TS_ASSERT( AnalysisDataService::Instance().doesExist("pdfE") );
- ITableWorkspace_sptr Etable = boost::dynamic_pointer_cast<ITableWorkspace>(
- API::AnalysisDataService::Instance().retrieve("pdfE"));
-
- TS_ASSERT(Etable);
- TS_ASSERT_EQUALS(Etable->columnCount(), 4);
- TS_ASSERT_EQUALS(Etable->rowCount(), n);
- TS_ASSERT_EQUALS(Etable->getColumn(0)->type(), "str");
- TS_ASSERT_EQUALS(Etable->getColumn(0)->name(), "Name");
- TS_ASSERT_EQUALS(Etable->getColumn(1)->type(), "double");
- TS_ASSERT_EQUALS(Etable->getColumn(1)->name(), "Value");
- TS_ASSERT_EQUALS(Etable->getColumn(2)->type(), "double");
- TS_ASSERT_EQUALS(Etable->getColumn(2)->name(), "Left's error");
- TS_ASSERT_EQUALS(Etable->getColumn(3)->type(), "double");
- TS_ASSERT_EQUALS(Etable->getColumn(3)->name(), "Rigth's error");
- TS_ASSERT(Etable->Double(0,1) == fun->getParameter("Height"));
- TS_ASSERT(Etable->Double(1,1) == fun->getParameter("Lifetime"));
+ TS_ASSERT( AnalysisDataService::Instance().doesExist("Parameters") );
+ ITableWorkspace_sptr Ptable = boost::dynamic_pointer_cast<ITableWorkspace>(
+ API::AnalysisDataService::Instance().retrieve("Parameters"));
+
+ TS_ASSERT(Ptable);
+ TS_ASSERT_EQUALS(Ptable->columnCount(), 4);
+ TS_ASSERT_EQUALS(Ptable->rowCount(), n);
+ TS_ASSERT_EQUALS(Ptable->getColumn(0)->type(), "str");
+ TS_ASSERT_EQUALS(Ptable->getColumn(0)->name(), "Name");
+ TS_ASSERT_EQUALS(Ptable->getColumn(1)->type(), "double");
+ TS_ASSERT_EQUALS(Ptable->getColumn(1)->name(), "Value");
+ TS_ASSERT_EQUALS(Ptable->getColumn(2)->type(), "double");
+ TS_ASSERT_EQUALS(Ptable->getColumn(2)->name(), "Left's error");
+ TS_ASSERT_EQUALS(Ptable->getColumn(3)->type(), "double");
+ TS_ASSERT_EQUALS(Ptable->getColumn(3)->name(), "Rigth's error");
+ TS_ASSERT(Ptable->Double(0,1) == fun->getParameter("Height"));
+ TS_ASSERT(Ptable->Double(1,1) == fun->getParameter("Lifetime"));
}
@@ -150,7 +149,7 @@ public:
fit.setProperty("WorkspaceIndex",0);
fit.setProperty("CreateOutput",true);
fit.setProperty("MaxIterations",10);
- fit.setProperty("Minimizer", "FABADA,ChainLength=5000,ConvergenceCriteria = 0.01, OutputWorkspaceConverged=conv");
+ fit.setProperty("Minimizer", "FABADA,ChainLength=5000,StepsBetweenValues=10,ConvergenceCriteria = 0.01");
TS_ASSERT_THROWS( fit.execute(), std::runtime_error );