plotasset.cpp

Go to the documentation of this file.

/*****************************************************************************
    NumeRe: Framework fuer Numerische Rechnungen
    Copyright (C) 2021  Erik Haenel et al.
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
 
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
******************************************************************************/
 
#include "plotasset.hpp"
 
mglData duplicatePoints(const mglData& _mData);
 
 
void PlotAsset::create(PlotType _t, size_t nDim, size_t nAxes, const std::vector<size_t>& samples, size_t nLayers)
{
    if (_t == PT_NONE || !nDim || samples.size() < nDim || samples.size() < nAxes || !nLayers)
    {
        type = PT_NONE;
        return;
    }
 
    for (size_t s = 0; s < samples.size(); s++)
    {
        if (!samples[s])
        {
            type = PT_NONE;
            return;
        }
    }
 
    type = _t;
    data.resize(nLayers);
 
    for (size_t l = 0; l < nLayers; l++)
    {
        data[l].first.Create(samples[XCOORD],
                             nDim >= 2 ? samples[YCOORD] : 1,
                             nDim == 3 ? samples[ZCOORD] : 1);
        data[l].second.Create(samples[XCOORD],
                              nDim >= 2 ? samples[YCOORD] : 1,
                              nDim == 3 ? samples[ZCOORD] : 1);
        data[l].first.Put(NAN);
        data[l].second.Put(NAN);
    }
 
    if (nAxes)
    {
        axes.resize(nAxes);
 
        for (size_t a = 0; a < nAxes; a++)
        {
            axes[a].Create(samples[a]);
        }
    }
}
 
 
void PlotAsset::writeData(const mu::value_type& val, size_t layer, size_t x, size_t y, size_t z)
{
    if (layer >= data.size()
        || x >= (size_t)data[layer].first.nx
        || y >= (size_t)data[layer].first.ny
        || z >= (size_t)data[layer].first.nz)
        throw SyntaxError(SyntaxError::PLOT_ERROR, "", "");
 
    data[layer].first.a[x + data[layer].first.nx*y + data[layer].first.nx*data[layer].first.ny*z] = mu::isinf(val) ? NAN : val.real();
    data[layer].second.a[x + data[layer].second.nx*y + data[layer].second.nx*data[layer].second.ny*z] = mu::isinf(val) ? NAN : val.imag();
}
 
 
void PlotAsset::writeAxis(double val, size_t pos, PlotCoords c)
{
    if (c >= axes.size())
        return;
 
    if (pos >= (size_t)axes[c].nx)
        throw SyntaxError(SyntaxError::PLOT_ERROR, "", "");
 
    axes[c].a[pos] = isinf(val) ? NAN : val;
}
 
 
void PlotAsset::duplicatePoints()
{
    for (size_t l = 0; l < data.size(); l++)
    {
        data[l].first = ::duplicatePoints(data[l].first);
        data[l].second = ::duplicatePoints(data[l].second);
    }
 
    for (size_t n = 0; n < axes.size(); n++)
    {
        axes[n] = ::duplicatePoints(axes[n]);
    }
}
 
 
void PlotAsset::removeNegativeValues(PlotCoords c)
{
    if (c < axes.size())
    {
        for (int n = 0; n < axes[c].GetNN(); n++)
        {
            axes[c].a[n] = axes[c].a[n] >= 0.0 ? axes[c].a[n] : NAN;
        }
    }
    else if (c - axes.size() < data.size())
    {
        for (int i = 0; i < data[c-axes.size()].first.GetNN(); i++)
        {
            data[c-axes.size()].first.a[i] = data[c-axes.size()].first.a[i] >= 0.0
                                             ? data[c-axes.size()].first.a[i] : NAN;
            data[c-axes.size()].second.a[i] = data[c-axes.size()].second.a[i] >= 0.0
                                              ? data[c-axes.size()].second.a[i] : NAN;
        }
    }
}
 
 
Interval PlotAsset::getAxisInterval(PlotCoords c) const
{
    if (c >= axes.size())
        return Interval();
 
    return Interval(axes[c].Minimal(), axes[c].Maximal());
}
 
 
IntervalSet PlotAsset::getDataIntervals(size_t layer) const
{
    IntervalSet ivl;
 
    if (layer >= data.size())
    {
        // Always return an "valid" intervalset
        ivl.intervals.resize(4);
        return ivl;
    }
 
    ivl.intervals.push_back(Interval(data[layer].first.Minimal(), data[layer].first.Maximal()));
    ivl.intervals.push_back(Interval(data[layer].second.Minimal(), data[layer].second.Maximal()));
    ivl.intervals.push_back(ivl[0].combine(ivl[1]));
    mglData normData = norm(layer);
    ivl.intervals.push_back(Interval(normData.Minimal(), normData.Maximal()));
    ivl.setNames({"Re", "Im", "ReIm", "AbsReIm"});
 
    return ivl;
}
 
 
mglData PlotAsset::norm(size_t layer) const
{
    if (layer >= data.size())
        return mglData();
 
    mglData absData = data[layer].first * data[layer].first + data[layer].second*data[layer].second;
 
    for (int i = 0; i < absData.GetNN(); i++)
        absData.a[i] = sqrt(absData.a[i]);
 
    return absData;
}
 
 
mglData PlotAsset::arg(size_t layer) const
{
    if (layer >= data.size())
        return mglData();
 
    mglData phaseData(data[layer].first.nx, data[layer].first.ny, data[layer].first.nz);
 
    for (int i = 0; i < phaseData.GetNN(); i++)
        phaseData.a[i] = atan2(data[layer].second.a[i], data[layer].first.a[i]);
 
    return phaseData;
}
 
 
bool PlotAsset::isComplex(size_t layer) const
{
    if (layer >= data.size())
        return false;
 
    return data[layer].second.Minimal() != 0.0 || data[layer].second.Maximal() != 0.0;
}
 
 
void PlotAsset::applyModulus(PlotCoords c, double mod)
{
    if (c < axes.size())
    {
        for (int n = 0; n < axes[c].nx; n++)
        {
            axes[c].a[n] = ::fmod(axes[c].a[n], mod);
 
            if (axes[c].a[n] < 0)
                axes[c].a[n] += mod;
        }
    }
    else if (c - axes.size() < data.size())
    {
        for (int n = 0; n < data[c-axes.size()].first.GetNN(); n++)
        {
            data[c-axes.size()].first.a[n] = ::fmod(data[c-axes.size()].first.a[n], mod);
            data[c-axes.size()].second.a[n] = ::fmod(data[c-axes.size()].second.a[n], mod);
 
            if (data[c-axes.size()].first.a[n] < 0)
                data[c-axes.size()].first.a[n] += mod;
 
            if (data[c-axes.size()].second.a[n] < 0)
                data[c-axes.size()].second.a[n] += mod;
        }
    }
}
 
 
mglData PlotAsset::vectorsToMatrix() const
{
    int maxdim = data[0].first.nx;
 
    // Find the largest vector
    for (size_t val = 1; val < data.size(); val++)
    {
        if (maxdim < data[val].first.nx)
            maxdim = data[val].first.nx;
    }
 
    // Create the target array and copy the data
    mglData _mData(maxdim, data.size());
 
    for (size_t val = 0; val < data.size(); val++)
    {
        for (int col = 0; col < data[val].first.nx; col++)
        {
            _mData.a[col + val * maxdim] = data[val].first.a[col];
        }
    }
 
    return _mData;
}
 
 
IntervalSet PlotAsset::getWeightedRanges(size_t layer, double dLowerPercentage, double dUpperPercentage) const
{
    IntervalSet ivl;
    ivl.intervals.resize(2);
 
    if (layer >= data.size())
        return ivl;
 
    mglData cp(data[layer].first);
    size_t nLength = qSortDouble(cp.a, cp.GetNN());
    size_t lowerpos = rint(fabs(1.0-dLowerPercentage)/2.0*(nLength));
    size_t upperpos = rint((dUpperPercentage/2.0+0.5)*(nLength));
 
    ivl[0].reset(cp.a[lowerpos], cp.a[upperpos]);
 
    cp = data[layer].second;
    nLength = qSortDouble(cp.a, cp.GetNN());
    lowerpos = rint(fabs(1.0-dLowerPercentage)/2.0*(nLength));
    upperpos = rint((dUpperPercentage/2.0+0.5)*(nLength));
 
    ivl[1].reset(cp.a[lowerpos], cp.a[upperpos]);
 
    return ivl;
}
 
 
 
// PLOTASSETMANAGER
 
 
 
IntervalSet PlotAssetManager::getIntervalsOfType(PlotType t, int coord) const
{
    IntervalSet ivl;
    ivl.intervals.resize(DATIVLCOUNT);
 
    for (size_t i = 0; i < assets.size(); i++)
    {
        if (assets[i].type != t)
            continue;
 
        if (coord == ALLRANGES
            || (coord == ONLYLEFT && assets[i].boundAxes.find('l') != std::string::npos)
            || (coord == ONLYRIGHT && assets[i].boundAxes.find('r') != std::string::npos))
        {
            IntervalSet dataIntervals;
 
            for (size_t l = 0; l < assets[i].getLayers(); l++)
            {
                dataIntervals = assets[i].getDataIntervals(l);
 
                for (int n = 0; n < DATIVLCOUNT; n++)
                    ivl[n] = ivl[n].combine(dataIntervals[n]);
            }
        }
        else
        {
            if ((int)assets[i].getLayers() <= coord)
                continue;
 
            IntervalSet dataIntervals = assets[i].getDataIntervals(coord);
 
            for (int n = 0; n < DATIVLCOUNT; n++)
                ivl[n] = ivl[n].combine(dataIntervals[n]);
        }
    }
 
    return ivl;
}
 
 
IntervalSet PlotAssetManager::getAxisIntervalsOfType(PlotType t) const
{
    IntervalSet ivl;
    ivl.intervals.resize(3);
 
    for (size_t i = 0; i < assets.size(); i++)
    {
        if (assets[i].type != t)
            continue;
 
        for (size_t c = 0; c < std::min(3u, assets[i].getDim()); c++)
        {
            ivl[c] = ivl[c].combine(assets[i].getAxisInterval((PlotCoords)c));
        }
    }
 
    return ivl;
}
 
 
void PlotAssetManager::normalize(int t_animate)
{
    std::pair<double, double> maxnorm(0, 0);
 
    for (const PlotAsset& ass : assets)
    {
        for (size_t layer = 0; layer < ass.getLayers(); layer++)
        {
            IntervalSet ivl = ass.getDataIntervals(layer);
 
            if (maxnorm.first < std::max(fabs(ivl[REAL].min()), fabs(ivl[REAL].max())))
                maxnorm.first = std::max(fabs(ivl[REAL].min()), fabs(ivl[REAL].max()));
 
            if (maxnorm.second < std::max(fabs(ivl[IMAG].min()), fabs(ivl[IMAG].max())))
                maxnorm.second = std::max(fabs(ivl[IMAG].min()), fabs(ivl[IMAG].max()));
        }
    }
 
    if (!t_animate)
        m_maxnorm = maxnorm;
 
    for (PlotAsset& ass : assets)
    {
        for (size_t layer = 0; layer < ass.getLayers(); layer++)
        {
            if (maxnorm.first)
                ass.data[layer].first /= maxnorm.first * maxnorm.first / m_maxnorm.first;
 
            if (maxnorm.second)
                ass.data[layer].second /= maxnorm.second * maxnorm.second / m_maxnorm.second;
        }
    }
}
 
 
IntervalSet PlotAssetManager::getDataIntervals(int coord) const
{
    return getIntervalsOfType(PT_DATA, coord);
}
 
 
IntervalSet PlotAssetManager::getFunctionIntervals(int coord) const
{
    return getIntervalsOfType(PT_FUNCTION, coord);
}
 
 
IntervalSet PlotAssetManager::getWeightedFunctionIntervals(int coord, double dLowerPercentage, double dUpperPercentage) const
{
    IntervalSet ivl;
    ivl.intervals.resize(2);
 
    for (size_t i = 0; i < assets.size(); i++)
    {
        if (assets[i].type != PT_FUNCTION)
            continue;
 
        if (coord == ALLRANGES
            || (coord == ONLYLEFT && assets[i].boundAxes.find('l') != std::string::npos)
            || (coord == ONLYRIGHT && assets[i].boundAxes.find('r') != std::string::npos))
        {
            IntervalSet dataIntervals;
 
            for (size_t l = 0; l < assets[i].getLayers(); l++)
            {
                dataIntervals = assets[i].getWeightedRanges(l, dLowerPercentage, dUpperPercentage);
                ivl[REAL] = ivl[REAL].combine(dataIntervals[REAL]);
                ivl[IMAG] = ivl[IMAG].combine(dataIntervals[IMAG]);
            }
        }
        else
        {
            if ((int)assets[i].getLayers() <= coord)
                continue;
 
            IntervalSet dataIntervals = assets[i].getWeightedRanges(coord, dLowerPercentage, dUpperPercentage);
            ivl[REAL] = ivl[REAL].combine(dataIntervals[REAL]);
            ivl[IMAG] = ivl[IMAG].combine(dataIntervals[IMAG]);
        }
    }
 
    return ivl;
}
 
 
IntervalSet PlotAssetManager::getDataAxes() const
{
    return getAxisIntervalsOfType(PT_DATA);
}
 
 
IntervalSet PlotAssetManager::getFunctionAxes() const
{
    return getAxisIntervalsOfType(PT_FUNCTION);
}
 
 
void PlotAssetManager::weightedRange(int coord, Interval& ivl) const
{
    if (log(ivl.range()) > 5)
    {
        const double dPercentage = 0.975;
        const double dSinglePercentageValue = 0.99;
        double dSinglePercentageUse;
 
        if (coord == ALLRANGES)
            dSinglePercentageUse = dSinglePercentageValue;
        else
            dSinglePercentageUse = dPercentage;
 
        if (log(fabs(ivl.min())) <= 1)
            ivl = getWeightedFunctionIntervals(coord, 1.0, dSinglePercentageUse)[0];
        else if (log(fabs(ivl.max())) <= 1)
            ivl = getWeightedFunctionIntervals(coord, dSinglePercentageUse, 1.0)[0];
        else
            ivl = getWeightedFunctionIntervals(coord, dPercentage, dPercentage)[0];
    }
}
 
 
bool PlotAssetManager::hasDataPlots() const
{
    for (const PlotAsset& ass : assets)
    {
        if (ass.type == PT_DATA)
            return true;
    }
 
    return false;
}
 
 
void PlotAssetManager::applyCoordSys(CoordinateSystem coords, size_t every)
{
    // Do not apply cartesian coordinates
    if (coords == CARTESIAN)
        return;
 
    for (size_t i = 0; i < assets.size(); i += every)
    {
        // Special case for 1D plots, as there are only two
        // possible valid combinations
        if (assets[i].getDim() == 1)
        {
            switch (coords)
            {
                case POLAR_PZ:
                case SPHERICAL_PT:
                    assets[i].applyModulus(XCOORD, 2.0*M_PI);
                    assets[i].removeNegativeValues(YCOORD);
                    break;
                case POLAR_RP:
                case POLAR_RZ:
                case SPHERICAL_RP:
                case SPHERICAL_RT:
                    assets[i].applyModulus(YCOORD, 2.0*M_PI);
                    assets[i].removeNegativeValues(XCOORD);
                    break;
                case CARTESIAN: break;
            }
        }
        else
        {
            switch (coords)
            {
                case POLAR_PZ:
                    assets[i].applyModulus(XCOORD, 2.0*M_PI);
                    assets[i].removeNegativeValues(ZCOORD);
                    break;
                case POLAR_RP:
                    assets[i].applyModulus(YCOORD, 2.0*M_PI);
                    assets[i].removeNegativeValues(XCOORD);
                    break;
                case POLAR_RZ:
                    assets[i].applyModulus(ZCOORD, 2.0*M_PI);
                    assets[i].removeNegativeValues(XCOORD);
                    break;
                case SPHERICAL_PT:
                    assets[i].applyModulus(XCOORD, 2.0*M_PI);
                    assets[i].applyModulus(YCOORD, 1.00001*M_PI);
                    assets[i].removeNegativeValues(ZCOORD);
                    break;
                case SPHERICAL_RP:
                    assets[i].applyModulus(YCOORD, 2.0*M_PI);
                    assets[i].applyModulus(ZCOORD, 1.00001*M_PI);
                    assets[i].removeNegativeValues(XCOORD);
                    break;
                case SPHERICAL_RT:
                    assets[i].applyModulus(ZCOORD, 2.0*M_PI);
                    assets[i].applyModulus(YCOORD, 1.00001*M_PI);
                    assets[i].removeNegativeValues(XCOORD);
                    break;
                case CARTESIAN: break;
            }
        }
    }
}