SimulationToPython.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Sim/Export/SimulationToPython.cpp
//! @brief     Implements class SimulationToPython.
//!
//! @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/Export/SimulationToPython.h"
#include "Base/Axis/IAxis.h"
#include "Base/Py/PyFmt.h"
#include "Base/Util/Algorithms.h"
#include "Base/Util/Assert.h"
#include "Device/Beam/Beam.h"
#include "Device/Beam/FootprintGauss.h"
#include "Device/Beam/FootprintSquare.h"
#include "Device/Detector/RectangularDetector.h"
#include "Device/Detector/SphericalDetector.h"
#include "Device/Mask/DetectorMask.h"
#include "Device/Resolution/ConvolutionDetectorResolution.h"
#include "Device/Resolution/ResolutionFunction2DGaussian.h"
#include "Param/Distrib/Distributions.h"
#include "Param/Distrib/RangedDistributions.h"
#include "Param/Node/NodeUtils.h"
#include "Resample/Options/SimulationOptions.h"
#include "Sim/Background/ConstantBackground.h"
#include "Sim/Background/PoissonBackground.h"
#include "Sim/Export/PyFmt2.h"
#include "Sim/Export/PyFmtLimits.h"
#include "Sim/Export/SampleToPython.h"
#include "Sim/Scan/AlphaScan.h"
#include "Sim/Scan/QzScan.h"
#include "Sim/Scan/ScanResolution.h"
#include "Sim/Simulation/includeSimulations.h"
#include <iomanip>
 
using Py::Fmt::indent;
 
namespace {
 
//! Returns a function that converts a coordinate to a Python code snippet with appropriate unit
std::function<std::string(double)> printFunc(const IDetector& detector)
{
    if (detector.defaultCoords() == Coords::MM)
        return Py::Fmt::printDouble;
    if (detector.defaultCoords() == Coords::RADIANS)
        return Py::Fmt::printDegrees;
    throw std::runtime_error("SimulationToPython::defineMasks() -> Error. Unknown detector units.");
}
 
bool isQuadraticDetector(const IDetector& det)
{
    ASSERT(det.rank() == 2);
    if (det.axis(0).size() != det.axis(1).size())
        return false;
    if (std::abs(det.axis(0).span() - det.axis(1).span())
        > 1e-12 * (det.axis(0).span() + det.axis(1).span()))
        return false;
    return true;
}
 
//! Returns true if it is (0, -1, 0) vector
bool isDefaultDirection(const R3 direction)
{
    return BaseUtils::algo::almostEqual(direction.x(), 0.0)
           && BaseUtils::algo::almostEqual(direction.y(), -1.0)
           && BaseUtils::algo::almostEqual(direction.z(), 0.0);
}
 
std::string defineFootprintFactor(const IFootprintFactor& foot)
{
    std::ostringstream result;
    result << indent() << "footprint = ba." << foot.name();
    result << "(" << Py::Fmt::printDouble(foot.widthRatio()) << ")\n";
    return result.str();
}
 
std::string defineScanResolution(const ScanResolution& scan)
{
    std::ostringstream result;
    result << Py::Fmt2::printRangedDistribution(*scan.distribution()) << "\n"
           << indent() << "resolution = "
           << "ba." << scan.name() << "(distribution, " << Py::Fmt::printDouble(scan.delta())
           << ")\n";
    return result.str();
}
 
std::string defineAlphaScan(const AlphaScan& scan)
{
    std::ostringstream result;
    result << "\n"
           << indent() << "# Define specular scan:\n"
           << indent() << "axis = " << Py::Fmt2::printAxis(scan.coordinateAxis(), "rad") << "\n"
           << indent() << "scan = "
           << "ba.AlphaScan(" << Py::Fmt::printDouble(scan.wavelength()) << ", axis)\n";
 
    if (scan.footprintFactor()) {
        result << defineFootprintFactor(*scan.footprintFactor());
        result << indent() << "scan.setFootprintFactor(footprint)\n";
    }
    if (const auto* r = scan.angleResolution(); r && r->distribution()) {
        result << defineScanResolution(*r) << "\n";
        result << indent() << "scan.setAngleResolution(resolution)\n";
    }
    if (const auto* r = scan.wavelengthResolution(); r && r->distribution()) {
        result << defineScanResolution(*r) << "\n";
        result << indent() << "scan.setWavelengthResolution(resolution)\n";
    }
    return result.str();
}
 
std::string defineQzScan(const QzScan& scan)
{
    std::ostringstream result;
    const std::string axis_def = indent() + "axis = ";
    result << axis_def << Py::Fmt2::printAxis(scan.coordinateAxis(), "") << "\n";
    // TODO correct unit would be 1/nm
 
    result << indent() << "scan = ba.QzScan(axis)";
    if (scan.resolution()) {
        result << "\n";
        result << defineScanResolution(*scan.resolution()) << "\n";
        result << indent() << "scan.setQResolution(resolution)";
    }
    return result.str();
}
 
std::string definePolarizationAnalyzer(const PolFilter& analyzer, const std::string parent)
{
    std::ostringstream result;
    R3 analyzer_direction = analyzer.analyzerDirection();
    double analyzer_efficiency = analyzer.analyzerEfficiency();
    double analyzer_total_transmission = analyzer.totalTransmission();
 
    if (analyzer_direction.mag() > 0.0) {
        std::string direction_name = "analyzer_direction";
        result << indent() << direction_name << " = R3("
               << Py::Fmt::printDouble(analyzer_direction.x()) << ", "
               << Py::Fmt::printDouble(analyzer_direction.y()) << ", "
               << Py::Fmt::printDouble(analyzer_direction.z()) << ")\n";
        result << indent() << parent << ".setAnalyzer(" << direction_name << ", "
               << Py::Fmt::printDouble(analyzer_efficiency) << ", "
               << Py::Fmt::printDouble(analyzer_total_transmission) << ")\n";
    }
    return result.str();
}
 
std::string defineScan(const ISpecularScan* scan)
{
    if (const auto* s = dynamic_cast<const AlphaScan*>(scan); s)
        return defineAlphaScan(*s);
    if (const auto* s = dynamic_cast<const QzScan*>(scan); s)
        return defineQzScan(*s);
    ASSERT(0);
}
 
std::string defineDetector(const ISimulation2D* simulation)
{
    const IDetector& detector = simulation->detector();
    if (detector.rank() != 2)
        throw std::runtime_error("defineDetector: "
                                 "detector must be two-dimensional for GISAS");
    std::ostringstream result;
    result << std::setprecision(12);
 
    if (const auto* const det = dynamic_cast<const SphericalDetector*>(&detector)) {
        ASSERT(det->rank() == 2);
        result << indent() << "detector = ba.SphericalDetector(";
        if (isQuadraticDetector(*det)) {
            result << det->axis(0).size() << ", " << Py::Fmt::printDegrees(det->axis(0).span())
                   << ", " << Py::Fmt::printDegrees(det->axis(0).center()) << ", "
                   << Py::Fmt::printDegrees(det->axis(1).center());
        } else {
            result << det->axis(0).size() << ", " << Py::Fmt::printDegrees(det->axis(0).min())
                   << ", " << Py::Fmt::printDegrees(det->axis(0).max()) << ", "
                   << det->axis(1).size() << ", " << Py::Fmt::printDegrees(det->axis(1).min())
                   << ", " << Py::Fmt::printDegrees(det->axis(1).max());
        }
        result << ")\n";
    } else if (const auto* const det = dynamic_cast<const RectangularDetector*>(&detector)) {
        result << "\n";
        result << indent() << "detector = ba.RectangularDetector(" << det->xSize() << ", "
               << Py::Fmt::printDouble(det->width()) << ", " << det->ySize() << ", "
               << Py::Fmt::printDouble(det->height()) << ")\n";
        if (det->getDetectorArrangment() == RectangularDetector::GENERIC) {
            result << indent() << "detector.setDetectorPosition("
                   << Py::Fmt::printKvector(det->getNormalVector()) << ", "
                   << Py::Fmt::printDouble(det->getU0()) << ", "
                   << Py::Fmt::printDouble(det->getV0());
            if (!isDefaultDirection(det->getDirectionVector()))
                result << ", " << Py::Fmt::printKvector(det->getDirectionVector());
            result << ")\n";
        } else if (det->getDetectorArrangment() == RectangularDetector::PERPENDICULAR_TO_SAMPLE) {
            result << indent() << "detector.setPerpendicularToSampleX("
                   << Py::Fmt::printDouble(det->getDistance()) << ", "
                   << Py::Fmt::printDouble(det->getU0()) << ", "
                   << Py::Fmt::printDouble(det->getV0()) << ")\n";
        } else if (det->getDetectorArrangment()
                   == RectangularDetector::PERPENDICULAR_TO_DIRECT_BEAM) {
            result << indent() << "detector.setPerpendicularToDirectBeam("
                   << Py::Fmt::printDouble(det->getDistance()) << ", "
                   << Py::Fmt::printDouble(det->getU0()) << ", "
                   << Py::Fmt::printDouble(det->getV0()) << ")\n";
        } else if (det->getDetectorArrangment()
                   == RectangularDetector::PERPENDICULAR_TO_REFLECTED_BEAM) {
            result << indent() << "detector.setPerpendicularToReflectedBeam("
                   << Py::Fmt::printDouble(det->getDistance()) << ", "
                   << Py::Fmt::printDouble(det->getU0()) << ", "
                   << Py::Fmt::printDouble(det->getV0()) << ")\n";
        } else if (det->getDetectorArrangment()
                   == RectangularDetector::PERPENDICULAR_TO_REFLECTED_BEAM_DPOS) {
            result << indent() << "detector.setPerpendicularToReflectedBeam("
                   << Py::Fmt::printDouble(det->getDistance()) << ")\n";
            result << indent() << "detector.setDirectBeamPosition("
                   << Py::Fmt::printDouble(det->getDirectBeamU0()) << ", "
                   << Py::Fmt::printDouble(det->getDirectBeamV0()) << ")\n";
        } else
            throw std::runtime_error("defineDetector() -> Error. Unknown alignment.");
    } else
        throw std::runtime_error("defineDetector() -> Error. Unknown detector");
    if (detector.hasExplicitRegionOfInterest()) {
        const auto xBounds = detector.regionOfInterestBounds(0);
        const auto yBounds = detector.regionOfInterestBounds(1);
        result << indent() << "detector.setRegionOfInterest(" << printFunc(detector)(xBounds.first)
               << ", " << printFunc(detector)(yBounds.first) << ", "
               << printFunc(detector)(xBounds.second) << ", " << printFunc(detector)(yBounds.second)
               << ")\n";
    }
    result << definePolarizationAnalyzer(detector.analyzer(), "detector");
    return result.str();
}
 
std::string defineDetectorResolutionFunction(const IDetector& detector)
{
    std::ostringstream result;
 
    if (const IDetectorResolution* resfunc = detector.detectorResolution()) {
        if (const auto* convfunc = dynamic_cast<const ConvolutionDetectorResolution*>(resfunc)) {
            if (const auto* resfunc = dynamic_cast<const ResolutionFunction2DGaussian*>(
                    convfunc->getResolutionFunction2D())) {
                result << indent() << "simulation.detector().setResolutionFunction(";
                result << "ba.ResolutionFunction2DGaussian(";
                result << printFunc(detector)(resfunc->sigmaX()) << ", ";
                result << printFunc(detector)(resfunc->sigmaY()) << "))\n";
            } else
                throw std::runtime_error("defineDetectorResolutionFunction() -> Error. "
                                         "Unknown detector resolution function");
        } else
            throw std::runtime_error("defineDetectorResolutionFunction() -> Error. "
                                     "Not a ConvolutionDetectorResolution function");
    }
    return result.str();
}
 
std::string defineBeamPolarization(const Beam& beam)
{
    std::ostringstream result;
    auto bloch_vector = beam.polVector();
    if (bloch_vector.mag() > 0.0) {
        std::string beam_polMatrices = "beam_polMatrices";
        result << indent() << beam_polMatrices << " = R3(" << Py::Fmt::printDouble(bloch_vector.x())
               << ", " << Py::Fmt::printDouble(bloch_vector.y()) << ", "
               << Py::Fmt::printDouble(bloch_vector.z()) << ")\n";
        result << indent() << "beam.setPolarization(" << beam_polMatrices << ")\n";
    }
    return result.str();
}
 
std::string defineBeamIntensity(const Beam& beam)
{
    std::ostringstream result;
    double beam_intensity = beam.intensity();
    if (beam_intensity > 0.0)
        result << indent() << "simulation.beam().setIntensity("
               << Py::Fmt::printScientificDouble(beam_intensity) << ")\n";
    return result.str();
}
 
std::string defineGISASBeam(const ScatteringSimulation& simulation)
{
    std::ostringstream result;
    const Beam& beam = simulation.beam();
 
    result << indent() << "beam = ba.Beam(" << Py::Fmt::printDouble(beam.intensity()) << ", "
           << Py::Fmt::printNm(beam.wavelength()) << ", ba.Direction("
           << Py::Fmt::printDegrees(beam.direction().alpha()) << ", "
           << Py::Fmt::printDegrees(beam.direction().phi()) << "))\n";
 
    result << defineBeamPolarization(beam);
 
    return result.str();
}
 
std::string defineOffspecBeam(const OffspecSimulation& simulation)
{
    std::ostringstream result;
    const Beam& beam = simulation.beam();
 
    const std::string axidef = indent() + "alpha_i_axis = ";
    result << axidef << Py::Fmt2::printAxis(simulation.beamAxis(), "rad") << "\n";
 
    result << indent() << "simulation.setBeamParameters(" << Py::Fmt::printNm(beam.wavelength())
           << ", "
           << "alpha_i_axis, " << Py::Fmt::printDegrees(beam.direction().phi()) << ")\n";
 
    result << defineBeamPolarization(beam);
    result << defineBeamIntensity(beam);
    return result.str();
}
 
std::string defineSpecularScan(const SpecularSimulation& simulation)
{
    std::ostringstream result;
    const ISpecularScan* scan = simulation.scan();
    if (!scan)
        throw std::runtime_error("Error defineSpecularScan: passed simulation "
                                 "does not contain any scan");
    result << defineScan(scan) << "\n";
    const PolFilter* analyzer = scan->analyzer();
    if (analyzer)
        result << definePolarizationAnalyzer(*analyzer, "scan");
    return result.str();
}
 
std::string defineParameterDistributions(const ISimulation* simulation)
{
    std::ostringstream result;
    const std::vector<ParameterDistribution>& distributions = simulation->getDistributions();
    if (distributions.empty())
        return "";
    for (size_t i = 0; i < distributions.size(); ++i) {
        const std::string mainParUnits =
            simulation->unitOfParameter(distributions[i].whichParameter());
 
        const std::string distr = "distr_" + std::to_string(i + 1);
        result << indent() << distr << " = "
               << "ba." << distributions[i].getDistribution()->pythonConstructor(mainParUnits)
               << "\n";
 
        result << indent() << "simulation.addParameterDistribution(ba."
               << distributions[i].whichParameterAsPyEnum() << ", " << distr << ", "
               << distributions[i].nDraws() << ", "
               << Py::Fmt::printDouble(distributions[i].sigmaFactor())
               << Py::Fmt::printRealLimitsArg(distributions[i].getLimits(), mainParUnits) << ")\n";
    }
    return result.str();
}
 
std::string defineMasks(const IDetector& detector)
{
    std::ostringstream result;
    result << std::setprecision(12);
 
    const DetectorMask* detectorMask = detector.detectorMask();
    if (detectorMask && detectorMask->hasMasks()) {
        result << "\n";
        for (size_t i_mask = 0; i_mask < detectorMask->numberOfMasks(); ++i_mask) {
            const MaskPattern* pat = detectorMask->patternAt(i_mask);
            IShape2D* shape = pat->shape;
            bool mask_value = pat->doMask;
            result << Py::Fmt2::representShape2D(indent(), shape, mask_value, printFunc(detector));
        }
        result << "\n";
    }
    return result.str();
}
 
std::string defineSimulationOptions(const ISimulation* simulation)
{
    std::ostringstream result;
    result << std::setprecision(12);
 
    const SimulationOptions& options = simulation->options();
    if (options.getHardwareConcurrency() != options.getNumberOfThreads())
        result << indent() << "simulation.options().setNumberOfThreads("
               << options.getNumberOfThreads() << ")\n";
    if (options.isIntegrate())
        result << indent() << "simulation.options().setMonteCarloIntegration(True, "
               << options.getMcPoints() << ")\n";
    if (options.useAvgMaterials())
        result << indent() << "simulation.options().setUseAvgMaterials(True)\n";
    if (options.includeSpecular())
        result << indent() << "simulation.options().setIncludeSpecular(True)\n";
    return result.str();
}
 
std::string defineBackground(const ISimulation* simulation)
{
    std::ostringstream result;
 
    const auto* bg = simulation->background();
    if (const auto* constant_bg = dynamic_cast<const ConstantBackground*>(bg)) {
        if (constant_bg->backgroundValue() > 0.0) {
            result << indent() << "background = ba.ConstantBackground("
                   << Py::Fmt::printScientificDouble(constant_bg->backgroundValue()) << ")\n";
            result << indent() << "simulation.setBackground(background)\n";
        }
    } else if (dynamic_cast<const PoissonBackground*>(bg)) {
        result << indent() << "background = ba.PoissonBackground()\n";
        result << indent() << "simulation.setBackground(background)\n";
    }
    return result.str();
}
 
std::string defineScatteringSimulation(const ScatteringSimulation* simulation)
{
    std::ostringstream result;
    result << defineGISASBeam(*simulation);
    result << defineDetector(simulation);
    result << indent() << "simulation = ba.ScatteringSimulation(beam, sample, detector)\n";
    result << defineDetectorResolutionFunction(simulation->detector());
    result << defineParameterDistributions(simulation);
    result << defineMasks(simulation->detector());
    result << defineSimulationOptions(simulation);
    result << defineBackground(simulation);
    return result.str();
}
 
std::string defineOffspecSimulation(const OffspecSimulation* simulation)
{
    std::ostringstream result;
    result << indent() << "simulation = ba.OffspecSimulation()\n";
    result << defineDetector(simulation);
    result << defineDetectorResolutionFunction(simulation->detector());
    result << defineOffspecBeam(*simulation);
    result << defineParameterDistributions(simulation);
    result << defineMasks(simulation->detector());
    result << defineSimulationOptions(simulation);
    result << defineBackground(simulation);
    return result.str();
}
 
std::string defineSpecularSimulation(const SpecularSimulation* simulation)
{
    std::ostringstream result;
    result << defineSpecularScan(*simulation);
    result << indent() << "simulation = ba.SpecularSimulation(scan, sample)\n";
    result << defineParameterDistributions(simulation);
    result << defineSimulationOptions(simulation);
    result << defineBackground(simulation);
    // result << defineBeamIntensity(simulation->beam());
    result << "\n";
    return result.str();
}
 
std::string defineSimulate(const ISimulation* simulation)
{
    std::ostringstream result;
    result << "def get_simulation(sample):\n";
    if (const auto* gisas = dynamic_cast<const ScatteringSimulation*>(simulation))
        result << defineScatteringSimulation(gisas);
    else if (const auto* offspec = dynamic_cast<const OffspecSimulation*>(simulation))
        result << defineOffspecSimulation(offspec);
    else if (const auto* spec = dynamic_cast<const SpecularSimulation*>(simulation))
        result << defineSpecularSimulation(spec);
    else
        ASSERT(0);
    result << "    return simulation\n\n\n";
 
    return result.str();
}
 
std::string simulationCode(const ISimulation& simulation)
{
    if (simulation.sample() == nullptr)
        throw std::runtime_error("Cannot export: Simulation has no sample");
    std::string code =
        SampleToPython().sampleCode(*simulation.sample()) + defineSimulate(&simulation);
    return "import bornagain as ba\n" + Py::Fmt::printImportedSymbols(code) + "\n\n" + code;
}
 
} // namespace
 
//  ************************************************************************************************
//  class SimulationToPython
//  ************************************************************************************************
 
std::string SimulationToPython::simulationPlotCode(const ISimulation& simulation)
{
    return simulationCode(simulation)
           + "if __name__ == '__main__':\n"
             "    from bornagain import ba_plot as bp\n"
             "    bp.parse_args()\n"
             "    sample = get_sample()\n"
             "    simulation = get_simulation(sample)\n"
             "    result = simulation.simulate()\n"
             "    bp.plot_simulation_result(result)\n";
}
 
std::string SimulationToPython::simulationSaveCode(const ISimulation& simulation,
                                                   const std::string& fname)
{
    return simulationCode(simulation)
           + "if __name__ == '__main__':\n"
             "    sample = get_sample()\n"
             "    simulation = get_simulation(sample)\n"
             "    result = simulation.simulate()\n"
             "    ba.IOFactory.writeSimulationResult(result, \""
           + fname + "\")\n";
}