Added class to change parameter values according to distributions;

refactored other classes to accomodate this change

Core/Algorithms/inc/Beam.h

@@ -39,10 +39,7 @@ public:
     virtual ~Beam() {}
 
     //! Get the value of the wavevector
-    cvector_t getCentralK() const { return m_central_k; }
-
-    //! Sets the value of the incoming wavevector
-    void setCentralK(const cvector_t& k_i);
+    cvector_t getCentralK() const;
 
     //! Sets the value of the incoming wavevector in terms of wavelength
     //! and incoming angles
@@ -84,7 +81,7 @@ private:
     //! Initialize polarization (for constructors)
     void initPolarization();
 
-    cvector_t m_central_k;  //!< incoming wavevector
+    double m_lambda, m_alpha, m_phi; //!< wavelength and angles of beam
     double m_intensity;     //!< beam intensity (neutrons/sec)
 #ifndef GCCXML_SKIP_THIS
     Eigen::Matrix2cd m_polarization; //!< polarization density matrix

Core/Algorithms/inc/DistributionWeighter.h

+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      Algorithms/inc/DistributionWeighter.h
+//! @brief     Defines class DistributionWeighter.
+//!
+//! @homepage  http://apps.jcns.fz-juelich.de/BornAgain
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2013
+//! @authors   Scientific Computing Group at MLZ Garching
+//! @authors   C. Durniak, G. Pospelov, W. Van Herck, J. Wuttke
+//
+// ************************************************************************** //
+
+#ifndef DISTRIBUTIONWEIGHTER_H_
+#define DISTRIBUTIONWEIGHTER_H_
+
+#include "IParameterized.h"
+#include "ParameterDistribution.h"
+
+
+class IDistribution1D;
+
+//! @class DistributionWeighter
+//! @ingroup algorithms_internal
+//! @brief Provides the functionality to average over parameter distributions
+//! with weights
+class DistributionWeighter : public IParameterized
+{
+public:
+	DistributionWeighter();
+	~DistributionWeighter();
+
+	//! add a samples parameter distribution
+	void addParameterDistribution(const std::string &param_name,
+			const IDistribution1D &distribution, size_t nbr_samples,
+			double sigma_factor=0.0);
+
+	//! get the total number of parameter value combinations (product
+	//! of the individual sizes of each parameter distribution
+	size_t getTotalNumberOfSamples() const;
+
+	//! set the parameter values of the simulation object to a specific
+	//! combination of values, determined by the index (which must be smaller
+	//! than the total number of combinations) and returns the weight
+	//! associated with this combination of parameter values
+	double setParameterValues(ParameterPool *p_parameter_pool, size_t index);
+private:
+	size_t m_nbr_combinations;
+	std::vector<ParameterDistribution> m_distributions;
+	std::vector<std::vector<ParameterSample> > m_cached_samples;
+};
+
+
+#endif /* DISTRIBUTIONWEIGHTER_H_ */

Core/Algorithms/src/Beam.cpp

@@ -18,15 +18,20 @@
 #include <Eigen/LU>
 
 
-Beam::Beam() : m_central_k(), m_intensity(1.0)
+Beam::Beam()
+: m_lambda(1.0)
+, m_alpha(0.0)
+, m_phi(0.0)
+, m_intensity(1.0)
 {
     setName("Beam");
     init_parameters();
     initPolarization();
 }
 
-Beam::Beam(const Beam& other) : IParameterized(), m_central_k(other.m_central_k)
-  , m_intensity(other.m_intensity), m_polarization(other.m_polarization)
+Beam::Beam(const Beam& other) : IParameterized(), m_lambda(other.m_lambda),
+		m_alpha(other.m_alpha), m_phi(other.m_phi),
+		m_intensity(other.m_intensity), m_polarization(other.m_polarization)
 {
     setName(other.getName());
     init_parameters();
@@ -42,17 +47,19 @@ Beam& Beam::operator=(const Beam& other)
     return *this;
 }
 
-void Beam::setCentralK(const cvector_t& k_i)
+cvector_t Beam::getCentralK() const
 {
-    m_central_k = k_i;
+    cvector_t k;
+    k.setLambdaAlphaPhi(m_lambda, m_alpha, m_phi);
+    return k;
 }
 
 void Beam::setCentralK(double lambda, double alpha_i, double phi_i)
 {
-    cvector_t k_i;
     if (alpha_i >0) alpha_i = - alpha_i;
-    k_i.setLambdaAlphaPhi(lambda, alpha_i, phi_i);
-    m_central_k = k_i;
+    m_lambda = lambda;
+    m_alpha = alpha_i;
+    m_phi = phi_i;
 }
 
 void Beam::setPolarization(const Eigen::Matrix2cd &polarization)
@@ -91,11 +98,16 @@ void Beam::init_parameters()
 {
     clearParameterPool();
     registerParameter("intensity", &m_intensity);
+    registerParameter("wavelength", &m_lambda);
+    registerParameter("alpha", &m_alpha);
+    registerParameter("phi", &m_phi);
 }
 
 void Beam::swapContent(Beam& other)
 {
-    std::swap(this->m_central_k, other.m_central_k);
+    std::swap(this->m_lambda, other.m_lambda);
+    std::swap(this->m_alpha, other.m_alpha);
+    std::swap(this->m_phi, other.m_phi);
     std::swap(this->m_intensity, other.m_intensity);
     std::swap(this->m_polarization, other.m_polarization);
 }
@@ -106,3 +118,4 @@ void Beam::initPolarization()
     m_polarization(0,0) = 0.5;
     m_polarization(1,1) = 0.5;
 }
+

Core/Algorithms/src/DistributionWeighter.cpp

+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      Algorithms/src/DistributionWeighter.cpp
+//! @brief     Implements class DistributionWeighter.
+//!
+//! @homepage  http://apps.jcns.fz-juelich.de/BornAgain
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2013
+//! @authors   Scientific Computing Group at MLZ Garching
+//! @authors   C. Durniak, G. Pospelov, W. Van Herck, J. Wuttke
+//
+// ************************************************************************** //
+
+#include "DistributionWeighter.h"
+
+DistributionWeighter::DistributionWeighter()
+: m_nbr_combinations(1)
+{
+}
+
+DistributionWeighter::~DistributionWeighter()
+{
+}
+
+void DistributionWeighter::addParameterDistribution(
+		const std::string &param_name, const IDistribution1D &distribution,
+		size_t nbr_samples, double sigma_factor)
+{
+	if (nbr_samples > 0) {
+		ParameterDistribution par_distr(param_name);
+		par_distr.setDistribution(distribution, nbr_samples, sigma_factor);
+		m_distributions.push_back(par_distr);
+		m_nbr_combinations *= nbr_samples;
+		m_cached_samples.push_back(par_distr.generateSamples());
+	}
+}
+
+size_t DistributionWeighter::getTotalNumberOfSamples() const
+{
+	return m_nbr_combinations;
+}
+
+double DistributionWeighter::setParameterValues(ParameterPool *p_parameter_pool,
+		size_t index)
+{
+	if (index >= m_nbr_combinations) {
+		throw Exceptions::RuntimeErrorException(
+				"DistributionWeighter::setParameterValues: "
+				"index must be smaller than the total number of parameter "
+				"combinations");
+	}
+	size_t n_distr = m_distributions.size();
+	if (n_distr == 0) return 0.0;
+	double weight = 1.0;
+	for (size_t param_index=n_distr-1; ; --param_index) {
+		size_t remainder = index % m_distributions[param_index].getNbrSamples();
+		index /= m_distributions[param_index].getNbrSamples();
+		int changed = p_parameter_pool->setMatchedParametersValue(
+				m_distributions[param_index].getParameterName(),
+				m_cached_samples[param_index][remainder].value);
+		if (changed != 1) {
+			throw Exceptions::RuntimeErrorException(
+					"DistributionWeighter::setParameterValues: "
+					" parameter name matches nothing or more than "
+					"one parameter");
+		}
+		weight *= m_cached_samples[param_index][remainder].weight;
+		if (param_index==0) break;
+	}
+	return weight;
+}

Core/Tools/inc/IParameterized.h

@@ -39,10 +39,12 @@ public:
     //! Returns pointer to the parameter pool.
     const ParameterPool* getParameterPool() const { return& m_parameters; }
 
-    //! Creates new parameter pool, with all local parameter and parameters of children
+    //! Creates new parameter pool, with all local parameters
+    //! and parameters of children
     virtual ParameterPool* createParameterTree() const;
 
-    //! Adds parameters from local pool to external pool and call recursion over direct children.
+    //! Adds parameters from local pool to external pool and
+    //! call recursion over direct children.
     virtual std::string addParametersToExternalPool(
         std::string path,
         ParameterPool* external_pool,

Core/Tools/inc/ParameterDistribution.h

@@ -28,14 +28,23 @@ class IDistribution1D;
 class ParameterDistribution : public IParameterized
 {
 public:
-	ParameterDistribution(const std::string& name);
+	ParameterDistribution(const std::string &name);
+	ParameterDistribution(const ParameterDistribution &other);
 	~ParameterDistribution();
 
+	//! Overload assignment operator
+	ParameterDistribution& operator=(const ParameterDistribution &other);
+
 	//! set the distribution type, number of samples and
 	//! the range of sample values
 	void setDistribution(const IDistribution1D &distribution,
 			size_t nbr_samples, double sigma_factor=0.0);
 
+	//! get the parameter's name
+	std::string getParameterName() const {
+		return m_name;
+	}
+
 	//! get number of samples for this distribution
 	size_t getNbrSamples() const {
 		return m_nbr_samples;

Core/Tools/src/ParameterDistribution.cpp

@@ -27,10 +27,31 @@ ParameterDistribution::ParameterDistribution(const std::string& name)
 	init_parameters();
 }
 
+ParameterDistribution::ParameterDistribution(
+		const ParameterDistribution& other)
+: m_name(other.m_name)
+, m_nbr_samples(other.m_nbr_samples)
+, m_sigma_factor(other.m_sigma_factor)
+{
+	mP_distribution.reset(other.mP_distribution->clone());
+}
+
 ParameterDistribution::~ParameterDistribution()
 {
 }
 
+ParameterDistribution& ParameterDistribution::operator=(
+		const ParameterDistribution& other)
+{
+	if (this != &other) {
+		m_name = other.m_name;
+		m_nbr_samples = other.m_nbr_samples;
+		m_sigma_factor = other.m_sigma_factor;
+		mP_distribution.reset(other.mP_distribution->clone());
+	}
+	return *this;
+}
+
 void ParameterDistribution::setDistribution(const IDistribution1D& distribution,
 		size_t nbr_samples, double sigma_factor)
 {

Tests/UnitTests/TestCore/BeamTest.h

@@ -30,12 +30,16 @@ BeamTest::~BeamTest()
 
 TEST_F(BeamTest, BeamInitialState)
 {
-    EXPECT_EQ(complex_t(0.0,0.0), emptyBeam.getCentralK()[0]);
+    EXPECT_DOUBLE_EQ(2.0*M_PI, emptyBeam.getCentralK()[0].real());
+    EXPECT_DOUBLE_EQ(0.0, emptyBeam.getCentralK()[0].imag());
     EXPECT_EQ(complex_t(0.0,0.0), emptyBeam.getCentralK()[1]);
     EXPECT_EQ(complex_t(0.0,0.0), emptyBeam.getCentralK()[2]);
     EXPECT_EQ(double(1.0), emptyBeam.getIntensity());
-    EXPECT_EQ(size_t(1), emptyBeam.getParameterPool()->size());
+    EXPECT_EQ(size_t(4), emptyBeam.getParameterPool()->size());
     EXPECT_EQ(double(1.0), emptyBeam.getParameterPool()->getParameter("intensity").getValue() );
+    EXPECT_EQ(double(1.0), emptyBeam.getParameterPool()->getParameter("wavelength").getValue() );
+    EXPECT_EQ(double(0.0), emptyBeam.getParameterPool()->getParameter("alpha").getValue() );
+    EXPECT_EQ(double(0.0), emptyBeam.getParameterPool()->getParameter("phi").getValue() );
     EXPECT_EQ(complex_t(0.5,0.0), emptyBeam.getPolarization()(0,0));
     EXPECT_EQ(complex_t(0.5,0.0), emptyBeam.getPolarization()(1,1));
     EXPECT_TRUE(emptyBeam.checkPolarization(emptyBeam.getPolarization()));
@@ -61,7 +65,7 @@ TEST_F(BeamTest, BeamAssignment)
     EXPECT_NEAR(2.85664, assignedBeam.getCentralK()[1].real(), 0.00001);
     EXPECT_NEAR(-5.28712, assignedBeam.getCentralK()[2].real(), 0.00001);
     EXPECT_EQ(double(2.0), assignedBeam.getIntensity());
-    EXPECT_EQ(size_t(1), assignedBeam.getParameterPool()->size());
+    EXPECT_EQ(size_t(4), assignedBeam.getParameterPool()->size());
     EXPECT_EQ(double(2.0), assignedBeam.getParameterPool()->getParameter("intensity").getValue() );
     EXPECT_EQ(complex_t(0.6,0.0), assignedBeam.getPolarization()(0,0));
     EXPECT_EQ(complex_t(0.4,0.0), assignedBeam.getPolarization()(1,1));