IHistogram.cpp

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Device/Histo/IHistogram.cpp
//! @brief     Implements class IntensityDataObject.
//!
//! @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/Histo/IHistogram.h"
#include "Base/Axis/Frame.h"
#include "Base/Math/Numeric.h"
#include "Base/Util/Assert.h"
#include "Device/Data/Powerfield.h"
#include "Device/Histo/Histogram1D.h"
#include "Device/Histo/Histogram2D.h"
#include "Device/Histo/IOFactory.h"
#include <memory>

IHistogram::IHistogram(const IAxis& axis_x)
    : m_frame(new Frame({axis_x.clone()}))
    , m_data(new Powerfield<CumulativeValue>(axis_x))
{
}

IHistogram::IHistogram(const IAxis& axis_x, const IAxis& axis_y)
    : m_frame(new Frame({axis_x.clone(), axis_y.clone()}))
    , m_data(new Powerfield<CumulativeValue>(axis_x, axis_y))
{
}

IHistogram::IHistogram(Powerfield<CumulativeValue>* data)
    : m_frame(new Frame(data->cloned_axes()))
    , m_data(data)
{
}


IHistogram::IHistogram(const IHistogram& other)
{
    m_data.reset(other.m_data->clone());
}


size_t IHistogram::getTotalNumberOfBins() const
{
    return m_data->allocatedSize();
}

const IAxis& IHistogram::xAxis() const
{
    return m_frame->axis(0);
}

const IAxis& IHistogram::yAxis() const
{
    return m_frame->axis(1);
}

double IHistogram::xMin() const
{
    return xAxis().min();
}

double IHistogram::xMax() const
{
    return xAxis().max();
}

size_t IHistogram::getNbinsX() const
{
    return xAxis().size();
}

double IHistogram::yMin() const
{
    return yAxis().min();
}

double IHistogram::yMax() const
{
    return yAxis().max();
}

size_t IHistogram::getNbinsY() const
{
    return yAxis().size();
}

size_t IHistogram::getGlobalBin(size_t binx, size_t biny) const
{
    std::vector<unsigned> axes_indices;
    axes_indices.push_back(static_cast<unsigned>(binx));
    if (rank() == 2)
        axes_indices.push_back(static_cast<unsigned>(biny));
    return m_data->toGlobalIndex(axes_indices);
}

size_t IHistogram::findGlobalBin(double x, double y) const
{
    std::vector<double> coordinates;
    coordinates.push_back(x);
    if (rank() == 2)
        coordinates.push_back(y);
    return m_data->findGlobalIndex(coordinates);
}

size_t IHistogram::xAxisIndex(size_t i) const
{
    return m_frame->projectedIndex(i, 0);
}

size_t IHistogram::yAxisIndex(size_t i) const
{
    return m_frame->projectedIndex(i, 1);
}

double IHistogram::xAxisValue(size_t i)
{
    return m_frame->projectedCoord(i, 0);
}

double IHistogram::yAxisValue(size_t i)
{
    return m_frame->projectedCoord(i, 1);
}

double IHistogram::binContent(size_t i) const
{
    return (*m_data)[i].getContent();
}

double IHistogram::binContent(size_t binx, size_t biny) const
{
    return binContent(getGlobalBin(binx, biny));
}

void IHistogram::setBinContent(size_t i, double value)
{
    (*m_data)[i].setContent(value);
}

void IHistogram::addBinContent(size_t i, double value)
{
    (*m_data)[i].add(value);
}

double IHistogram::binError(size_t i) const
{
    return (*m_data)[i].getRMS();
}

double IHistogram::binError(size_t binx, size_t biny) const
{
    return binError(getGlobalBin(binx, biny));
}

double IHistogram::binAverage(size_t i) const
{
    return (*m_data)[i].getAverage();
}

double IHistogram::binAverage(size_t binx, size_t biny) const
{
    return binAverage(getGlobalBin(binx, biny));
}

double IHistogram::getMaximum() const
{
    ASSERT(m_data->allocatedSize()>0);
    double ret = (*m_data)[0].getContent();
    for (size_t i=1; i < m_data->allocatedSize(); ++i)
        if ((*m_data)[i].getContent() > ret)
            ret = (*m_data)[i].getContent();
    return ret;
}

double IHistogram::getMinimum() const
{
    ASSERT(m_data->allocatedSize()>0);
    double ret = (*m_data)[0].getContent();
    for (size_t i=1; i < m_data->allocatedSize(); ++i)
        if ((*m_data)[i].getContent() < ret)
            ret = (*m_data)[i].getContent();
    return ret;
}

void IHistogram::scale(double value)
{
    for (size_t index = 0; index < getTotalNumberOfBins(); ++index)
        (*m_data)[index].setContent(value * (*m_data)[index].getContent());
}

double IHistogram::integral() const
{
    double result(0.0);
    for (size_t index = 0; index < getTotalNumberOfBins(); ++index)
        result += (*m_data)[index].getContent();
    return result;
}

#ifdef BORNAGAIN_PYTHON
PyObject* IHistogram::array(DataType dataType) const
{
    const std::unique_ptr<Powerfield<double>> data(createPowerfield(dataType));
    return data->getArray();
}
#endif // BORNAGAIN_PYTHON

void IHistogram::reset()
{
    m_data->setAllTo(CumulativeValue());
}

IHistogram* IHistogram::createHistogram(const Powerfield<double>& source)
{
    if (source.rank() == 1)
        return new Histogram1D(source);
    if (source.rank() == 2)
        return new Histogram2D(source);
    std::ostringstream message;
    message << "IHistogram::createHistogram(const Powerfield<double>& source) -> Error. ";
    message << "The rank of source " << source.rank() << " ";
    message << "is not suitable for creation neither 1-dim nor 2-dim histograms.";
    throw std::runtime_error(message.str());
}

std::vector<double> IHistogram::meanVector() const
{
    std::vector<double> result(getTotalNumberOfBins());
    for (size_t i = 0; i < result.size(); ++i)
        result[i] = binContent(i);
    return result;
}

void IHistogram::init_from_data(const Powerfield<double>& source)
{
    if (rank() != source.rank()) {
        std::ostringstream message;
        message << "IHistogram::IHistogram(const Powerfield<double>& data) -> Error. ";
        message << "The dimension of this histogram " << rank() << " ";
        message << "is differ from the dimension of source " << m_data->rank() << std::endl;
        throw std::runtime_error(message.str());
    }

    m_data.reset(new Powerfield<CumulativeValue>(source.cloned_axes()));
    for (size_t i = 0; i < source.allocatedSize(); ++i)
        (*m_data)[i].add(source[i]);
}

//! Returns data of requested type for globalbin number
double IHistogram::binData(size_t i, IHistogram::DataType dataType) const
{
    if (dataType == DataType::INTEGRAL)
        return binContent(i);
    if (dataType == DataType::AVERAGE)
        return binAverage(i);
    if (dataType == DataType::STANDARD_ERROR)
        return binError(i);
    throw std::runtime_error("IHistogram::binData() -> Error. Unknown data type.");
}

//! Returns vector of values of requested DataType
std::vector<double> IHistogram::getDataVector(IHistogram::DataType dataType) const
{
    std::vector<double> result;
    result.resize(getTotalNumberOfBins(), 0.0);
    for (size_t index = 0; index < getTotalNumberOfBins(); ++index)
        result[index] = binData(index, dataType);
    return result;
}

//! Copy content (but not the axes) from other histogram. Dimensions should be the same.
void IHistogram::copyContentFrom(const IHistogram& other)
{
    if (!hasSameSizes(other))
        throw std::runtime_error(
            "IHistogram::copyContentFrom() -> Error. Can't copy the data of different shape.");
    reset();
    for (size_t i = 0; i < getTotalNumberOfBins(); ++i)
        (*m_data)[i] = (*other.m_data)[i];
}

//! creates new Powerfield with histogram's shape and put there values corresponding to DataType
Powerfield<double>* IHistogram::createPowerfield(IHistogram::DataType dataType) const
{
    auto* result = new Powerfield<double>(m_data->cloned_axes());
    for (size_t i = 0; i < getTotalNumberOfBins(); ++i)
        (*result)[i] = binData(i, dataType);
    return result;
}

bool IHistogram::hasSameShape(const IHistogram& other) const
{
    return m_data->hasSameShape(*other.m_data);
}

bool IHistogram::hasSameSizes(const IHistogram& other) const
{
    return m_data->hasSameSizes(*other.m_data);
}

const IHistogram& IHistogram::operator+=(const IHistogram& right)
{
    if (!hasSameSizes(right))
        throw std::runtime_error(
            "IHistogram::operator+=() -> Error. Histograms have different dimension");
    for (size_t i = 0; i < getTotalNumberOfBins(); ++i)
        addBinContent(i, right.binContent(i));
    return *this;
}

IHistogram* IHistogram::relativeDifferenceHistogram(const IHistogram& rhs) const
{
    if (!hasSameSizes(rhs))
        throw std::runtime_error("IHistogram::relativeDifferenceHistogram() -> Error. "
                                 "Histograms have different dimensions");

    IHistogram* result = this->clone();
    result->reset();

    for (size_t i = 0; i < getTotalNumberOfBins(); ++i) {
        double diff = Numeric::GetRelativeDifference(binContent(i), rhs.binContent(i));
        result->setBinContent(i, diff);
    }
    return result;
}

void IHistogram::save(const std::string& filename) const
{
    IOFactory::writeIntensityData(*this, filename);
}

void IHistogram::load(const std::string& filename)
{
    const std::unique_ptr<IHistogram> hist(IOFactory::readIntensityData(filename));
    copyContentFrom(*hist);
}