diff --git a/Core/Algorithms/inc/Beam.h b/Core/Algorithms/inc/Beam.h
index 8afec34ae2578a860831f2dd969f2ef057a2ab50..11a58878050f6863441b7c1d382ab5c4f0c567b7 100644
--- a/Core/Algorithms/inc/Beam.h
+++ b/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
diff --git a/Core/Algorithms/inc/DistributionWeighter.h b/Core/Algorithms/inc/DistributionWeighter.h
new file mode 100644
index 0000000000000000000000000000000000000000..cc70bacd9411fc2af494c9adfa7c96061a421279
--- /dev/null
+++ b/Core/Algorithms/inc/DistributionWeighter.h
@@ -0,0 +1,56 @@
+// ************************************************************************** //
+//
+// 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 ¶m_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_ */
diff --git a/Core/Algorithms/src/Beam.cpp b/Core/Algorithms/src/Beam.cpp
index e6c7e34197cb343995a0702c8cd1f6979d2b7d79..f549968a3bb581deb551a5474e7a05a84f538120 100644
--- a/Core/Algorithms/src/Beam.cpp
+++ b/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;
}
+
diff --git a/Core/Algorithms/src/DistributionWeighter.cpp b/Core/Algorithms/src/DistributionWeighter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..db9c6f9a10ca8127d1fc4226b442d89e71b10c86
--- /dev/null
+++ b/Core/Algorithms/src/DistributionWeighter.cpp
@@ -0,0 +1,73 @@
+// ************************************************************************** //
+//
+// 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 ¶m_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;
+}
diff --git a/Core/Tools/inc/IParameterized.h b/Core/Tools/inc/IParameterized.h
index 87b5974de4602bb47c29eda166983c54d8c11471..81172e9182a55c7e7ba8f38cc6fbdd109d931722 100644
--- a/Core/Tools/inc/IParameterized.h
+++ b/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,
diff --git a/Core/Tools/inc/ParameterDistribution.h b/Core/Tools/inc/ParameterDistribution.h
index 3130ba3d1642dc05e11a8978644ba70259fbccfe..3ed787ed6a6ac85d0a923aea64d5721ce52b061c 100644
--- a/Core/Tools/inc/ParameterDistribution.h
+++ b/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;
diff --git a/Core/Tools/src/ParameterDistribution.cpp b/Core/Tools/src/ParameterDistribution.cpp
index 74a0e64a045787b292e81dc452020ec6264eadb5..38a74cddbdc77b7fdb0a2f30785715aa68e6cedb 100644
--- a/Core/Tools/src/ParameterDistribution.cpp
+++ b/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)
{
diff --git a/Tests/UnitTests/TestCore/BeamTest.h b/Tests/UnitTests/TestCore/BeamTest.h
index 04bd0d8467dd48eff855d94d68815b9b9918f942..15b5423b7c20ed946b452052d261377c5eecf84f 100644
--- a/Tests/UnitTests/TestCore/BeamTest.h
+++ b/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));