ISimulation.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Sim/Simulation/ISimulation.cpp
//! @brief     Implements interface ISimulation.
//!
//! @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/Simulation/ISimulation.h"
#include "Base/Progress/ProgressHandler.h"
#include "Base/Util/Assert.h"
#include "Base/Util/StringUtils.h"
#include "Device/Histo/SimulationResult.h"
#include "Resample/Options/SimulationOptions.h"
#include "Resample/Processed/ReSample.h"
#include "Sample/Multilayer/MultiLayer.h"
#include "Sample/Multilayer/MultilayerUtils.h"
#include "Sim/Background/IBackground.h"
#include "Sim/Computation/IComputation.h"
#include <gsl/gsl_errno.h>
#include <iostream>
#include <thread>
 
namespace {
 
size_t indexStep(size_t total_size, size_t n_handlers)
{
    ASSERT(total_size > 0);
    ASSERT(n_handlers > 0);
    size_t result = total_size / n_handlers;
    return total_size % n_handlers ? ++result : result;
}
 
size_t startIndex(size_t n_handlers, size_t current_handler, size_t n_elements)
{
    const size_t handler_size = indexStep(n_elements, n_handlers);
    const size_t start_index = current_handler * handler_size;
    if (start_index >= n_elements)
        return n_elements;
    return start_index;
}
 
size_t batchSize(size_t n_handlers, size_t current_handler, size_t n_elements)
{
    const size_t handler_size = indexStep(n_elements, n_handlers);
    const size_t start_index = current_handler * handler_size;
    if (start_index >= n_elements)
        return 0;
    return std::min(handler_size, n_elements - start_index);
}
 
void runComputations(const std::vector<std::unique_ptr<IComputation>>& computations)
{
    ASSERT(!computations.empty());
 
    if (computations.size() == 1) { // Running computation in current thread
        const auto& computation = computations.front();
        computation->compute();
        if (computation->isCompleted())
            return;
        std::string message = computation->errorMessage();
        throw std::runtime_error("Error in runComputations: ISimulation has "
                                 "terminated unexpectedly with following error: "
                                 "message.\n"
                                 + message);
    }
 
    // Running computations in several threads.
    // The number of threads is equal to the number of computations.
 
    std::vector<std::unique_ptr<std::thread>> threads;
 
    // Run simulations in n threads.
    for (const auto& comp : computations)
        threads.emplace_back(new std::thread([&comp]() { comp->compute(); }));
 
    // Wait for threads to complete.
    for (auto& thread : threads)
        thread->join();
 
    // Check successful completion.
    std::vector<std::string> failure_messages;
    for (const auto& comp : computations)
        if (!comp->isCompleted())
            failure_messages.push_back(comp->errorMessage());
 
    if (failure_messages.empty())
        return;
    throw std::runtime_error("Error in runComputations: "
                             "At least one simulation thread has terminated unexpectedly.\n"
                             "Messages: "
                             + BaseUtils::String::join(failure_messages, " --- "));
}
 
} // namespace
 
 
//  ************************************************************************************************
//  class implementation
//  ************************************************************************************************
 
ISimulation::ISimulation(const MultiLayer& sample)
    : m_sample(sample.clone())
    , m_options(std::make_unique<SimulationOptions>())
    , m_progress(std::make_unique<ProgressHandler>())
{
}
 
ISimulation::~ISimulation() = default;
 
const SimulationOptions& ISimulation::options() const
{
    ASSERT(m_options);
    return *m_options;
}
 
SimulationOptions& ISimulation::options()
{
    ASSERT(m_options);
    return *m_options;
}
 
ProgressHandler& ISimulation::progress()
{
    ASSERT(m_progress);
    return *m_progress;
}
 
void ISimulation::subscribe(const std::function<bool(size_t)>& inform)
{
    ASSERT(m_progress);
    m_progress->subscribe(inform);
}
 
//! Initializes a progress monitor that prints to stdout.
void ISimulation::setTerminalProgressMonitor()
{
    m_progress->subscribe([](size_t percentage_done) -> bool {
        if (percentage_done < 100)
            std::cout << std::setprecision(2) << "\r... " << percentage_done << "%" << std::flush;
        else // wipe out
            std::cout << "\r... 100%\n";
        return true;
    });
}
 
//! Runs simulation with possible averaging over parameter distributions; returns result.
SimulationResult ISimulation::simulate()
{
    if (!SampleUtils::Multilayer::ContainsCompatibleMaterials(*m_sample))
        throw std::runtime_error("Sample contains incompatible material definitions");
    gsl_set_error_handler_off();
 
    prepareSimulation();
 
    const auto re_sample = reSample::make(*sample(), options(), force_polarized());
 
    const size_t total_size = numberOfElements();
    size_t param_combinations = m_distribution_handler.getTotalNumberOfSamples();
 
    m_progress->reset();
    m_progress->setExpectedNTicks(param_combinations * total_size);
 
    // restrict calculation to current batch
    const size_t n_batches = m_options->getNumberOfBatches();
    const size_t current_batch = m_options->getCurrentBatch();
 
    const size_t batch_start = startIndex(n_batches, current_batch, total_size);
    const size_t batch_size = batchSize(n_batches, current_batch, total_size);
    ASSERT(batch_size);
 
    if (param_combinations == 1)
        runSingleSimulation(re_sample, batch_start, batch_size, 1.);
    else {
        initDistributionHandler();
        for (size_t index = 0; index < param_combinations; ++index) {
            double weight = m_distribution_handler.setParameterValues(index);
            runSingleSimulation(re_sample, batch_start, batch_size, weight);
        }
    }
    moveDataFromCache();
    return pack_result();
}
 
const MultiLayer* ISimulation::sample() const
{
    return m_sample.get();
}
 
void ISimulation::setBackground(const IBackground& bg)
{
    m_background.reset(bg.clone());
}
 
std::vector<const INode*> ISimulation::nodeChildren() const
{
    std::vector<const INode*> result;
    if (m_sample)
        result << m_sample.get();
    return result;
}
 
std::string ISimulation::unitOfParameter(ParameterDistribution::WhichParameter which) const
{
    switch (which) {
    case ParameterDistribution::BeamWavelength:
        return "nm";
    case ParameterDistribution::BeamAzimuthalAngle:
    case ParameterDistribution::BeamInclinationAngle:
        return "rad";
    default:
        throw std::runtime_error("ISimulation::unitOfParameter(): missing implementation");
    }
}
 
void ISimulation::addParameterDistribution(ParameterDistribution::WhichParameter whichParameter,
                                           const IDistribution1D& distribution, size_t nbr_samples,
                                           double sigma_factor, const RealLimits& limits)
{
    ParameterDistribution par_distr(whichParameter, distribution, nbr_samples, sigma_factor,
                                    limits);
    addParameterDistribution(par_distr);
}
 
void ISimulation::addParameterDistribution(const ParameterDistribution& par_distr)
{
    validateParametrization(par_distr);
    m_distribution_handler.addParameterDistribution(par_distr);
}
 
const std::vector<ParameterDistribution>& ISimulation::getDistributions() const
{
    return m_distribution_handler.getDistributions();
}
 
 
//! Runs a single simulation with fixed parameter values.
//! If desired, the simulation is run in several threads.
void ISimulation::runSingleSimulation(const reSample& re_sample, size_t batch_start,
                                      size_t batch_size, double weight)
{
    initElementVector();
 
    const size_t n_threads = m_options->getNumberOfThreads();
    ASSERT(n_threads > 0);
 
    std::vector<std::unique_ptr<IComputation>> computations;
 
    for (size_t i_thread = 0; i_thread < n_threads; ++i_thread) {
        const size_t thread_start = batch_start + startIndex(n_threads, i_thread, batch_size);
        const size_t thread_size = batchSize(n_threads, i_thread, batch_size);
        if (thread_size == 0)
            break;
        computations.emplace_back(createComputation(re_sample, thread_start, thread_size));
    }
    runComputations(computations);
 
    normalize(batch_start, batch_size);
    addBackgroundIntensity(batch_start, batch_size);
    addDataToCache(weight);
}