SimDataPair.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Sim/Fitting/SimDataPair.cpp
//! @brief     Defines class SimDataPair.
//!
//! @homepage  http://www.bornagainproject.org
//! @license   GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2018
//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
//  ************************************************************************************************
 
#include "Sim/Fitting/SimDataPair.h"
#include "Base/Axis/Frame.h"
#include "Base/Axis/IAxis.h"
#include "Base/Math/Numeric.h"
#include "Base/Util/Assert.h"
#include "Device/Coord/CoordSystem1D.h"
#include "Device/Coord/CoordSystem2D.h"
#include "Device/Data/DataUtils.h"
#include "Device/Data/Datafield.h"
#include "Device/Detector/IDetector.h"
#include "Sim/Scan/ISpecularScan.h"
#include "Sim/Simulation/ISimulation2D.h"
#include <utility>
 
namespace {
 
[[noreturn]] void throwInitializationException(std::string method)
{
    std::stringstream ss;
    ss << "Error in SimDataPair::" << method << ": Trying access non-initialized data\n";
    throw std::runtime_error(ss.str());
}
 
std::unique_ptr<Datafield> initUserWeights(const Datafield& shape, double value)
{
    auto result = std::make_unique<Datafield>(shape.frame().cloned_axes());
    result->setAllTo(value);
    return result;
}
 
bool haveSameSizes(const IDetector& detector, const Datafield& data)
{
    if (data.rank() != detector.rank())
        return false;
 
    for (size_t i = 0; i < detector.rank(); ++i)
        if (data.axis(i).size() != detector.axis(i).size())
            return false;
 
    return true;
}
 
//! Convert user data to SimulationResult object for later drawing in various axes units.
//! User data will be cropped to the ROI defined in the simulation, amplitudes in areas
//! corresponding to the masked areas of the detector will be set to zero.
 
SimulationResult convertData(const ISimulation2D& simulation, const Datafield& data)
{
    const ICoordSystem* coordSystem = simulation.createCoordSystem();
    auto roi_data = std::make_unique<Datafield>(coordSystem->defaultAxes());
 
    if (roi_data->hasSameSizes(data)) {
        // data is already cropped to ROI
        simulation.detector().iterateOverNonMaskedPoints([&](IDetector::const_iterator it) {
            (*roi_data)[it.roiIndex()] = data[it.roiIndex()];
        });
    } else if (haveSameSizes(simulation.detector(), data)) {
        // experimental data has same shape as the detector, we have to copy the original
        // data to a smaller roi map
        simulation.detector().iterateOverNonMaskedPoints([&](IDetector::const_iterator it) {
            (*roi_data)[it.roiIndex()] = data[it.detectorIndex()];
        });
    } else
        throw std::runtime_error(
            "FitObject::init_dataset: Detector and experimental data have different shape");
 
    return SimulationResult(*roi_data, std::move(coordSystem));
}
 
} // namespace
 
//  ************************************************************************************************
//  class implementation
//  ************************************************************************************************
 
SimDataPair::SimDataPair(simulation_builder_t builder, const Datafield& raw_data,
                         std::unique_ptr<Datafield>&& raw_stdv, double user_weight)
    : m_simulation_builder(std::move(builder))
    , m_raw_data(raw_data.clone())
    , m_raw_uncertainties(std::move(raw_stdv))
{
    m_raw_user_weights = initUserWeights(*m_raw_data, user_weight);
    validate();
}
 
SimDataPair::SimDataPair(simulation_builder_t builder, const Datafield& raw_data,
                         std::unique_ptr<Datafield>&& raw_stdv,
                         std::unique_ptr<Datafield>&& user_weights)
    : m_simulation_builder(std::move(builder))
    , m_raw_data(raw_data.clone())
    , m_raw_uncertainties(std::move(raw_stdv))
    , m_raw_user_weights(std::move(user_weights))
{
    if (!m_raw_user_weights)
        m_raw_user_weights = initUserWeights(*m_raw_data, 1.0);
    validate();
}
 
SimDataPair::SimDataPair(SimDataPair&& other)
    : m_simulation_builder(std::move(other.m_simulation_builder))
    , m_sim_data(std::move(other.m_sim_data))
    , m_exp_data(std::move(other.m_exp_data))
    , m_uncertainties(std::move(other.m_uncertainties))
    , m_user_weights(std::move(other.m_user_weights))
    , m_raw_data(std::move(other.m_raw_data))
    , m_raw_uncertainties(std::move(other.m_raw_uncertainties))
    , m_raw_user_weights(std::move(other.m_raw_user_weights))
{
    validate();
}
 
SimDataPair::~SimDataPair() = default;
 
void SimDataPair::execSimulation(const mumufit::Parameters& params)
{
    std::unique_ptr<ISimulation> simulation = m_simulation_builder(params);
    ASSERT(simulation);
    m_sim_data = simulation->simulate();
 
    if (!m_exp_data.empty() && !m_uncertainties.empty() && !m_user_weights.empty())
        return;
 
    if (m_sim_data.empty())
        throwInitializationException("initResultArrays");
 
    auto* const sim2d = dynamic_cast<ISimulation2D*>(simulation.get());
    if (sim2d) {
        m_exp_data = convertData(*sim2d, *m_raw_data);
        m_user_weights = convertData(*sim2d, *m_raw_user_weights);
    } else {
        const ICoordSystem& converter = m_sim_data.converter();
        m_exp_data = SimulationResult(*m_raw_data, converter);
        m_user_weights = SimulationResult(*m_raw_user_weights, converter);
    }
 
    if (sim2d && containsUncertainties())
        m_uncertainties = convertData(*sim2d, *m_raw_uncertainties);
    else {
        const ICoordSystem& converter = m_sim_data.converter();
        auto dummy_array = std::make_unique<Datafield>(converter.defaultAxes());
        m_uncertainties = SimulationResult(*dummy_array, converter);
    }
}
 
bool SimDataPair::containsUncertainties() const
{
    return static_cast<bool>(m_raw_uncertainties);
}
 
SimulationResult SimDataPair::simulationResult() const
{
    if (m_sim_data.empty())
        throwInitializationException("simulationResult");
    return m_sim_data;
}
 
SimulationResult SimDataPair::experimentalData() const
{
    if (m_exp_data.empty())
        throwInitializationException("experimentalData");
    return m_exp_data;
}
 
SimulationResult SimDataPair::uncertainties() const
{
    if (m_uncertainties.empty())
        throwInitializationException("uncertainties");
    return m_uncertainties;
}
 
//! Returns the user uncertainties cut to the ROI area.
SimulationResult SimDataPair::userWeights() const
{
    if (m_user_weights.empty())
        throwInitializationException("userWeights");
    return m_user_weights;
}
 
//! Returns relative difference between simulation and experimental data.
 
SimulationResult SimDataPair::relativeDifference() const
{
    if (m_sim_data.empty() || m_exp_data.empty())
        throwInitializationException("relativeDifference");
 
    SimulationResult result = m_sim_data;
    for (size_t i = 0, size = result.size(); i < size; ++i)
        result[i] = Numeric::relativeDifference(result[i], m_exp_data[i]);
 
    return result;
}
 
SimulationResult SimDataPair::absoluteDifference() const
{
    if (m_sim_data.empty() || m_exp_data.empty())
        throwInitializationException("absoluteDifference");
 
    SimulationResult result = m_sim_data;
    for (size_t i = 0, size = result.size(); i < size; ++i)
        result[i] = Numeric::GetAbsoluteDifference(result[i], m_exp_data[i]);
 
    return result;
}
 
std::vector<double> SimDataPair::experimental_array() const
{
    if (m_exp_data.empty())
        throwInitializationException("experimental_array");
    return m_exp_data.datafield()->flatVector();
}
 
std::vector<double> SimDataPair::simulation_array() const
{
    if (m_sim_data.empty())
        throwInitializationException("simulation_array");
    return m_sim_data.datafield()->flatVector();
}
 
std::vector<double> SimDataPair::uncertainties_array() const
{
    if (m_uncertainties.empty())
        throwInitializationException("uncertainties_array");
    return m_uncertainties.datafield()->flatVector();
}
 
std::vector<double> SimDataPair::user_weights_array() const
{
    if (m_user_weights.empty())
        throwInitializationException("user_weights_array");
    return m_user_weights.datafield()->flatVector();
}
 
void SimDataPair::validate() const
{
    if (!m_simulation_builder)
        throw std::runtime_error("Error in SimDataPair: simulation builder is empty");
 
    if (!m_raw_data)
        throw std::runtime_error("Error in SimDataPair: passed experimental data array is empty");
 
    if (m_raw_uncertainties && !m_raw_uncertainties->hasSameShape(*m_raw_data))
        throw std::runtime_error(
            "Error in SimDataPair: experimental data and uncertainties have different shape.");
 
    if (!m_raw_user_weights || !m_raw_user_weights->hasSameShape(*m_raw_data))
        throw std::runtime_error(
            "Error in SimDataPair: user weights are not initialized or have invalid shape");
}