Datafield.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Device/Data/Datafield.cpp
//! @brief     Implements template specializations of class Datafield.cpp.
//!
//! @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 "Device/Data/Datafield.h"
#include "Base/Axis/FixedBinAxis.h"
#include "Base/Axis/Frame.h"
#include "Base/Py/PyCore.h"
#include "Base/Util/Assert.h"
#include "Base/Util/ThreadInfo.h"
#include <algorithm>
 
 
Datafield::Datafield(Frame* frame)
    : m_frame(frame)
    , m_vec(m_frame->size())
{
}
 
Datafield::Datafield(Frame* frame, const std::vector<double>& vector)
    : m_frame(frame)
    , m_vec(vector)
{
}
 
Datafield::Datafield(const std::vector<IAxis*>& axes, const std::vector<double>& vector)
    : Datafield(new Frame(axes), vector)
{
}
 
Datafield::Datafield(const std::vector<IAxis*>& axes)
    : Datafield(new Frame(axes))
{
}
 
Datafield::~Datafield() {}
 
Datafield* Datafield::clone() const
{
    return new Datafield(m_frame->cloned_axes(), m_vec);
}
 
void Datafield::setAt(size_t i, double val)
{
    m_vec[i] = val;
}
 
double Datafield::valAt(size_t i) const
{
    return m_vec[i];
}
 
void Datafield::setAllTo(const double& value)
{
    for (double& e : m_vec)
        e = value;
}
 
void Datafield::setVector(const std::vector<double>& vector)
{
    ASSERT(vector.size() == frame().size());
    m_vec = vector;
}
 
size_t Datafield::rank() const
{
    return m_frame->rank();
}
 
size_t Datafield::size() const
{
    ASSERT(frame().size() == m_vec.size());
    return frame().size();
}
 
const Frame& Datafield::frame() const
{
    return *m_frame;
}
 
const IAxis& Datafield::axis(size_t k) const
{
    return m_frame->axis(k);
}
 
const IAxis& Datafield::xAxis() const
{
    return m_frame->axis(0);
}
 
const IAxis& Datafield::yAxis() const
{
    return m_frame->axis(1);
}
 
 
//! Returns true if object have same dimensions
bool Datafield::hasSameSizes(const Datafield& other) const
{
    return frame().hasSameSizes(other.frame());
}
 
//! Returns true if object have same dimensions and shape of axis
bool Datafield::hasSameShape(const Datafield& other) const
{
    return frame() == other.frame();
}
 
std::vector<double> Datafield::flatVector() const
{
    std::vector<double> ret(size());
    for (size_t i = 0; i < size(); ++i)
        ret[i] = (*this)[i];
    return ret;
}
 
Datafield* Datafield::normalizedToMaximum() const
{
    double maxval = maxVal();
    ASSERT(maxval > 0);
    std::vector<double> out(frame().size());
    for (size_t i = 0; i < frame().size(); ++i)
        out[i] = m_vec[i] / maxval;
    return new Datafield(m_frame->cloned_axes(), out);
}
 
double Datafield::maxVal() const
{
    return *std::max_element(m_vec.begin(), m_vec.end());
}
 
double Datafield::minVal() const
{
    return *std::min_element(m_vec.begin(), m_vec.end());
}
 
Datafield* Datafield::crop(double xmin, double ymin, double xmax, double ymax) const
{
    const std::unique_ptr<IAxis> xaxis{xAxis().clone()};
    const std::unique_ptr<IAxis> yaxis{yAxis().clone()};
    xaxis->clip(xmin, xmax);
    yaxis->clip(ymin, ymax);
 
    std::vector<double> out(size());
    size_t iout = 0;
    for (size_t i = 0; i < size(); ++i) {
        double x = m_frame->projectedCoord(i, 0);
        double y = m_frame->projectedCoord(i, 1);
        if (xaxis->contains(x) && yaxis->contains(y))
            out[iout++] = m_vec[i];
    }
    return new Datafield(m_frame->cloned_axes(), out);
}
 
Datafield* Datafield::crop(double xmin, double xmax) const
{
    const std::unique_ptr<IAxis> xaxis{xAxis().clone()};
    xaxis->clip(xmin, xmax);
    std::vector<double> out(size());
    size_t iout = 0;
    for (size_t i = 0; i < size(); ++i) {
        const double x = m_frame->projectedCoord(i, 0);
        if (xaxis->contains(x))
            out[iout++] = m_vec[i];
    }
    return new Datafield(m_frame->cloned_axes(), out);
}
 
#ifdef BORNAGAIN_PYTHON
 
PyObject* Datafield::npArray() const
{
    std::vector<size_t> dimensions;
    for (size_t i = 0; i < rank(); i++)
        dimensions.push_back(axis(i).size());
 
    // for rot90 of 2-dim arrays to conform with numpy
    if (dimensions.size() == 2)
        std::swap(dimensions[0], dimensions[1]);
 
    // creating ndarray objects describing size of dimensions
    npy_int ndim_numpy = (int)dimensions.size();
    auto* ndimsizes_numpy = new npy_intp[dimensions.size()];
    for (size_t i = 0; i < dimensions.size(); i++)
        ndimsizes_numpy[i] = dimensions[i];
 
    // creating standalone numpy array
    PyObject* pyarray = PyArray_SimpleNew(ndim_numpy, ndimsizes_numpy, NPY_DOUBLE);
    delete[] ndimsizes_numpy;
    ASSERT(pyarray);
 
    // getting pointer to data buffer of numpy array
    double* array_buffer = (double*)PyArray_DATA((PyArrayObject*)pyarray);
 
    // filling numpy array with output_data
    if (rank() == 2) {
        for (size_t i = 0; i < size(); ++i) {
            std::vector<int> axes_indices = m_frame->allIndices(i);
            size_t offset =
                axes_indices[0] + axis(0).size() * (axis(1).size() - 1 - axes_indices[1]);
            array_buffer[offset] = (*this)[i];
        }
    } else {
        for (size_t i = 0; i < size(); ++i)
            *array_buffer++ = (*this)[i];
    }
 
    return pyarray;
}
 
#endif // BORNAGAIN_PYTHON
 
 
Datafield* Datafield::xProjection()
{
    return create_xProjection(0, static_cast<int>(xAxis().size()) - 1);
}
 
Datafield* Datafield::xProjection(double yvalue)
{
    int ybin_selected = static_cast<int>(yAxis().findClosestIndex(yvalue));
    return create_xProjection(ybin_selected, ybin_selected);
}
 
Datafield* Datafield::xProjection(double ylow, double yup)
{
    int ybinlow = static_cast<int>(yAxis().findClosestIndex(ylow));
    int ybinup = static_cast<int>(yAxis().findClosestIndex(yup));
    return create_xProjection(ybinlow, ybinup);
}
 
Datafield* Datafield::yProjection()
{
    return create_yProjection(0, static_cast<int>(xAxis().size()) - 1);
}
 
Datafield* Datafield::yProjection(double xvalue)
{
    int xbin_selected = static_cast<int>(xAxis().findClosestIndex(xvalue));
    return create_yProjection(xbin_selected, xbin_selected);
}
 
Datafield* Datafield::yProjection(double xlow, double xup)
{
    int xbinlow = static_cast<int>(xAxis().findClosestIndex(xlow));
    int xbinup = static_cast<int>(xAxis().findClosestIndex(xup));
    return create_yProjection(xbinlow, xbinup);
}
 
Datafield* Datafield::create_xProjection(int ybinlow, int ybinup)
{
    std::vector<double> out(xAxis().size());
    for (size_t i = 0; i < size(); ++i) {
        int ybin = static_cast<int>(m_frame->projectedIndex(i, 1));
        if (ybin >= ybinlow && ybin <= ybinup) {
            double x = m_frame->projectedCoord(i, 0);
            ASSERT(xAxis().contains(x));
            size_t iout = xAxis().findClosestIndex(x);
            out[iout] += valAt(i);
        }
    }
    return new Datafield({xAxis().clone()}, out);
}
 
Datafield* Datafield::create_yProjection(int xbinlow, int xbinup)
{
    std::vector<double> out(yAxis().size());
    for (size_t i = 0; i < size(); ++i) {
        int xbin = static_cast<int>(m_frame->projectedIndex(i, 0));
        if (xbin >= xbinlow && xbin <= xbinup) {
 
            // TODO: duplicates code from create_xProjection, move to Frame ?
 
            double y = m_frame->projectedCoord(i, 1);
            ASSERT(yAxis().contains(y));
            size_t iout = yAxis().findClosestIndex(y);
            out[iout] += valAt(i);
        }
    }
    return new Datafield({yAxis().clone()}, out);
}