XYZVertexTwoParticle.hpp

Go to the documentation of this file.

 
#pragma once
#include <istream>
#include "lib/ValueBundle.hpp"
#include "lib/Assert.hpp"
#include "FrgCommon.hpp"
 
template <int size> struct XYZVertexTwoParticleAccessBuffer
{
public:
    XYZVertexTwoParticleAccessBuffer() : siteExchange(false) {}
 
    int frequencyOffsets[size]; 
    float frequencyWeights[size]; 
    int signFlag[size]; 
    bool siteExchange; 
};
 
struct XYZVertexTwoParticle
{
public:
    enum struct FrequencyChannel
    {
        S, 
        T, 
        U, 
        All, 
        None 
    };  
 
    XYZVertexTwoParticle()
    {
        //store width in all memory dimensions
        _memoryStepLattice = FrgCommon::lattice().size;
        _memoryStepLatticeT = _memoryStepLattice * FrgCommon::frequency().size;
 
        sizeFrequency = FrgCommon::frequency().size * FrgCommon::frequency().size * (FrgCommon::frequency().size + 1) / 2;
        size = FrgCommon::lattice().size * sizeFrequency;
 
        //alloc and init memory
        _dataXX = new float[size];
        _dataYY = new float[size];
        _dataZZ = new float[size];
        _dataDD = new float[size];
        memset(_dataXX, 0, sizeof(float) * size);
        memset(_dataYY, 0, sizeof(float) * size);
        memset(_dataZZ, 0, sizeof(float) * size);
        memset(_dataDD, 0, sizeof(float) * size);
    }
 
    ~XYZVertexTwoParticle()
    {
        delete[] _dataXX;
        delete[] _dataYY;
        delete[] _dataZZ;
        delete[] _dataDD;
    }
 
    void expandIterator(int iterator, LatticeIterator &i1, float &s, float &t, float &u) const
    {
        ASSERT(iterator >= 0 && iterator < size);
        ASSERT(iterator >= 0 && iterator < size);
        ASSERT(&s != &t);
        ASSERT(&t != &u);
        ASSERT(&s != &u);
 
        int it = iterator;
        int su = it / _memoryStepLatticeT;
        it = it % _memoryStepLatticeT;
        t = FrgCommon::frequency()._data[it / _memoryStepLattice];
        i1 = FrgCommon::lattice().fromParametrization(it % _memoryStepLattice);
 
        for (int so = 0; so <= su; ++so)
        {
            for (int uo = 0; uo <= so; ++uo)
            {
                if (su == so * (so + 1) / 2 + uo)
                {
                    s = FrgCommon::frequency()._data[so];
                    u = FrgCommon::frequency()._data[uo];
                    return;
                }
            }
        }
    }
 
    void expandIterator(int iterator, float &s, float &t, float &u) const
    {
        ASSERT(iterator >= 0 && iterator < sizeFrequency);
        ASSERT(iterator >= 0 && iterator < size);
        ASSERT(&s != &t);
        ASSERT(&t != &u);
        ASSERT(&s != &u);
 
        int su = iterator / FrgCommon::frequency().size;
        t = FrgCommon::frequency()._data[iterator % FrgCommon::frequency().size];
 
        for (int so = 0; so <= su; ++so)
        {
            for (int uo = 0; uo <= so; ++uo)
            {
                if (su == so * (so + 1) / 2 + uo)
                {
                    s = FrgCommon::frequency()._data[so];
                    u = FrgCommon::frequency()._data[uo];
                    return;
                }
            }
        }
    }
 
    float &getValueRef(const int iterator, const SpinComponent symmetry) const
    {
        if (symmetry == SpinComponent::X) return _dataXX[iterator];
        if (symmetry == SpinComponent::Y) return _dataYY[iterator];
        if (symmetry == SpinComponent::Z) return _dataZZ[iterator];
        else return _dataDD[iterator];
    }
 
    float getValue(LatticeIterator i1, LatticeIterator i2, float s, float t, float u, SpinComponent symmetry, const XYZVertexTwoParticle::FrequencyChannel channel) const
    {
        //map to positive frequency sector
        if (s < 0 && u < 0)
        {
            s = -s;
            u = -u;
        }
        else
        {
            if (s < 0)
            {
                s = -s;
                std::swap(i1, i2);
            }
            else if (u < 0)
            {
                u = -u;
                std::swap(i1, i2);
            }
        }
        if (t < 0)
        {
            t = -t;
        }
 
        int siteOffset = FrgCommon::lattice().symmetryTransform(i1, i2, symmetry);
 
        if (channel == FrequencyChannel::S)
        {
            int exactS = FrgCommon::frequency().offset(s);
 
            int lowerT, upperT;
            float biasT;
            int lowerU, upperU;
            float biasU;
 
            FrgCommon::frequency().interpolateOffset(t, lowerT, upperT, biasT);
            FrgCommon::frequency().interpolateOffset(u, lowerU, upperU, biasU);
 
            return (1 - biasU) * (
                (1 - biasT) * (_directAccessMapFrequencyExchange(siteOffset, exactS, lowerT, lowerU, symmetry)) + biasT * (_directAccessMapFrequencyExchange(siteOffset, exactS, upperT, lowerU, symmetry))
                ) + biasU * (
                (1 - biasT) * (_directAccessMapFrequencyExchange(siteOffset, exactS, lowerT, upperU, symmetry)) + biasT * (_directAccessMapFrequencyExchange(siteOffset, exactS, upperT, upperU, symmetry))
                    );
        }
        else if (channel == FrequencyChannel::T)
        {
            int exactT = FrgCommon::frequency().offset(t);
 
            int lowerS, upperS;
            float biasS;
            int lowerU, upperU;
            float biasU;
 
            FrgCommon::frequency().interpolateOffset(s, lowerS, upperS, biasS);
            FrgCommon::frequency().interpolateOffset(u, lowerU, upperU, biasU);
 
            return (1 - biasU) * (
                (1 - biasS) * (_directAccessMapFrequencyExchange(siteOffset, lowerS, exactT, lowerU, symmetry)) + biasS * (_directAccessMapFrequencyExchange(siteOffset, upperS, exactT, lowerU, symmetry))
                ) + biasU * (
                (1 - biasS) * (_directAccessMapFrequencyExchange(siteOffset, lowerS, exactT, upperU, symmetry)) + biasS * (_directAccessMapFrequencyExchange(siteOffset, upperS, exactT, upperU, symmetry))
                    );
        }
        else if (channel == FrequencyChannel::U)
        {
            int exactU = FrgCommon::frequency().offset(u);
 
            int lowerS, upperS;
            float biasS;
            int lowerT, upperT;
            float biasT;
 
            FrgCommon::frequency().interpolateOffset(s, lowerS, upperS, biasS);
            FrgCommon::frequency().interpolateOffset(t, lowerT, upperT, biasT);
 
            return (1 - biasT) * (
                (1 - biasS) * (_directAccessMapFrequencyExchange(siteOffset, lowerS, lowerT, exactU, symmetry)) + biasS * (_directAccessMapFrequencyExchange(siteOffset, upperS, lowerT, exactU, symmetry))
                ) + biasT * (
                (1 - biasS) * (_directAccessMapFrequencyExchange(siteOffset, lowerS, upperT, exactU, symmetry)) + biasS * (_directAccessMapFrequencyExchange(siteOffset, upperS, upperT, exactU, symmetry))
                    );
        }
        else if (channel == FrequencyChannel::None)
        {
            int lowerS, upperS;
            float biasS;
            int lowerT, upperT;
            float biasT;
            int lowerU, upperU;
            float biasU;
 
            FrgCommon::frequency().interpolateOffset(s, lowerS, upperS, biasS);
            FrgCommon::frequency().interpolateOffset(t, lowerT, upperT, biasT);
            FrgCommon::frequency().interpolateOffset(u, lowerU, upperU, biasU);
 
            return
                (1 - biasU) * (
                (1 - biasT) * (
                    (1 - biasS) * (_directAccessMapFrequencyExchange(siteOffset, lowerS, lowerT, lowerU, symmetry)) + biasS * (_directAccessMapFrequencyExchange(siteOffset, upperS, lowerT, lowerU, symmetry))
                    ) + biasT * (
                    (1 - biasS) * (_directAccessMapFrequencyExchange(siteOffset, lowerS, upperT, lowerU, symmetry)) + biasS * (_directAccessMapFrequencyExchange(siteOffset, upperS, upperT, lowerU, symmetry))
                        )
                    ) + biasU * (
                    (1 - biasT) * (
                        (1 - biasS) * (_directAccessMapFrequencyExchange(siteOffset, lowerS, lowerT, upperU, symmetry)) + biasS * (_directAccessMapFrequencyExchange(siteOffset, upperS, lowerT, upperU, symmetry))
                        ) + biasT * (
                        (1 - biasS) * (_directAccessMapFrequencyExchange(siteOffset, lowerS, upperT, upperU, symmetry)) + biasS * (_directAccessMapFrequencyExchange(siteOffset, upperS, upperT, upperU, symmetry))
                            )
                        );
        }
        else if (channel == FrequencyChannel::All)
        {
            int exactS = FrgCommon::frequency().offset(s);
            int exactT = FrgCommon::frequency().offset(t);
            int exactU = FrgCommon::frequency().offset(u);
            return _directAccessMapFrequencyExchange(siteOffset, exactS, exactT, exactU, symmetry);
        }
        else
        {
            throw Exception(Exception::Type::ArgumentError, "Invalid interpolation channel");
            return 0.0f;
        }
    }
 
    template <int n> float getValue(const LatticeIterator i1, const LatticeIterator i2, SpinComponent symmetry, const XYZVertexTwoParticleAccessBuffer<n> &accessBuffer) const
    {
        int siteOffset = (accessBuffer.siteExchange) ? FrgCommon::lattice().symmetryTransform(i2, i1, symmetry) : FrgCommon::lattice().symmetryTransform(i1, i2, symmetry);
 
        float value = 0.0f;
        if (symmetry == SpinComponent::X)
        {
            for (int i = 0; i < n; ++i) value += accessBuffer.frequencyWeights[i] * _dataXX[accessBuffer.frequencyOffsets[i] + siteOffset];
        }
        else if (symmetry == SpinComponent::Y)
        {
            for (int i = 0; i < n; ++i) value += accessBuffer.frequencyWeights[i] * _dataYY[accessBuffer.frequencyOffsets[i] + siteOffset];
        }
        else if (symmetry == SpinComponent::Z)
        {
            for (int i = 0; i < n; ++i) value += accessBuffer.frequencyWeights[i] * _dataZZ[accessBuffer.frequencyOffsets[i] + siteOffset];
        }
        else
        {
            for (int i = 0; i < n; ++i) value += accessBuffer.signFlag[i] * accessBuffer.frequencyWeights[i] * _dataDD[accessBuffer.frequencyOffsets[i] + siteOffset];
        }
 
        return value;
    }
 
    template <int n> float getValueLocal(const SpinComponent symmetry, const XYZVertexTwoParticleAccessBuffer<n> &accessBuffer) const
    {
        float value = 0.0f;
        if (symmetry == SpinComponent::X)
        {
            for (int i = 0; i < n; ++i) value += accessBuffer.frequencyWeights[i] * _dataXX[accessBuffer.frequencyOffsets[i]];
        }
        else if (symmetry == SpinComponent::Y)
        {
            for (int i = 0; i < n; ++i) value += accessBuffer.frequencyWeights[i] * _dataYY[accessBuffer.frequencyOffsets[i]];
        }
        else if (symmetry == SpinComponent::Z)
        {
            for (int i = 0; i < n; ++i) value += accessBuffer.frequencyWeights[i] * _dataZZ[accessBuffer.frequencyOffsets[i]];
        }
        else
        {
            for (int i = 0; i < n; ++i) value += accessBuffer.signFlag[i] * accessBuffer.frequencyWeights[i] * _dataDD[accessBuffer.frequencyOffsets[i]];
        }
 
        return value;
    }
 
    template <int n> void getValueSuperbundle(const XYZVertexTwoParticleAccessBuffer<n> &accessBuffer, ValueSuperbundle<float, 4> &superbundle) const
    {
        superbundle.reset();
        const LatticeSiteDescriptor *sites = (accessBuffer.siteExchange) ? FrgCommon::lattice().getInvertedSites() : FrgCommon::lattice().getSites();
 
        float *base[4] = { _dataXX, _dataYY, _dataZZ, _dataDD };
 
        for (int i = 0; i < n; ++i)
        {
            float weight = accessBuffer.frequencyWeights[i];
            float signedWeight = accessBuffer.signFlag[i] * accessBuffer.frequencyWeights[i];
            int frequencyOffset = accessBuffer.frequencyOffsets[i];
            int size = FrgCommon::lattice().size;
 
            for (int j = 0; j < size; ++j)
            {
                superbundle.bundle(0)[j] += weight * base[static_cast<int>(sites[j].spinPermutation[0])][frequencyOffset + sites[j].rid];
                superbundle.bundle(1)[j] += weight * base[static_cast<int>(sites[j].spinPermutation[1])][frequencyOffset + sites[j].rid];
                superbundle.bundle(2)[j] += weight * base[static_cast<int>(sites[j].spinPermutation[2])][frequencyOffset + sites[j].rid];
                superbundle.bundle(3)[j] += signedWeight * base[3][frequencyOffset + sites[j].rid];
            }
        }
    }
 
    XYZVertexTwoParticleAccessBuffer<4> generateAccessBuffer(float s, float t, float u, const XYZVertexTwoParticle::FrequencyChannel channel) const
    {
        ASSERT(channel == FrequencyChannel::S || channel == FrequencyChannel::T || channel == FrequencyChannel::U);
 
        XYZVertexTwoParticleAccessBuffer<4> accessBuffer;
 
        //map to positive frequency sector
        if (s < 0 && u < 0)
        {
            s = -s;
            u = -u;
        }
        else
        {
            if (s < 0)
            {
                s = -s;
                accessBuffer.siteExchange = true;
            }
            else if (u < 0)
            {
                u = -u;
                accessBuffer.siteExchange = true;
            }
        }
        if (t < 0)
        {
            t = -t;
        }
 
        //interpolate frequency
        if (channel == FrequencyChannel::S)
        {
            int exactS = FrgCommon::frequency().offset(s);
            int lowerT, upperT;
            float biasT;
            int lowerU, upperU;
            float biasU;
            FrgCommon::frequency().interpolateOffset(t, lowerT, upperT, biasT);
            FrgCommon::frequency().interpolateOffset(u, lowerU, upperU, biasU);
 
            accessBuffer.frequencyWeights[0] = (1 - biasU) * (1 - biasT);
            accessBuffer.frequencyOffsets[0] = _generateAccessBufferOffset(exactS, lowerT, lowerU, accessBuffer.signFlag[0]);
            accessBuffer.frequencyWeights[1] = (1 - biasU) * biasT;
            accessBuffer.frequencyOffsets[1] = _generateAccessBufferOffset(exactS, upperT, lowerU, accessBuffer.signFlag[1]);
            accessBuffer.frequencyWeights[2] = biasU * (1 - biasT);
            accessBuffer.frequencyOffsets[2] = _generateAccessBufferOffset(exactS, lowerT, upperU, accessBuffer.signFlag[2]);
            accessBuffer.frequencyWeights[3] = biasU * biasT;
            accessBuffer.frequencyOffsets[3] = _generateAccessBufferOffset(exactS, upperT, upperU, accessBuffer.signFlag[3]);
        }
        else if (channel == FrequencyChannel::T)
        {
            int exactT = FrgCommon::frequency().offset(t);
            int lowerS, upperS;
            float biasS;
            int lowerU, upperU;
            float biasU;
            FrgCommon::frequency().interpolateOffset(s, lowerS, upperS, biasS);
            FrgCommon::frequency().interpolateOffset(u, lowerU, upperU, biasU);
 
            accessBuffer.frequencyWeights[0] = (1 - biasU) * (1 - biasS);
            accessBuffer.frequencyOffsets[0] = _generateAccessBufferOffset(lowerS, exactT, lowerU, accessBuffer.signFlag[0]);
            accessBuffer.frequencyWeights[1] = (1 - biasU) * biasS;
            accessBuffer.frequencyOffsets[1] = _generateAccessBufferOffset(upperS, exactT, lowerU, accessBuffer.signFlag[1]);
            accessBuffer.frequencyWeights[2] = biasU * (1 - biasS);
            accessBuffer.frequencyOffsets[2] = _generateAccessBufferOffset(lowerS, exactT, upperU, accessBuffer.signFlag[2]);
            accessBuffer.frequencyWeights[3] = biasU * biasS;
            accessBuffer.frequencyOffsets[3] = _generateAccessBufferOffset(upperS, exactT, upperU, accessBuffer.signFlag[3]);
 
        }
        else if (channel == FrequencyChannel::U)
        {
            int exactU = FrgCommon::frequency().offset(u);
            int lowerS, upperS;
            float biasS;
            int lowerT, upperT;
            float biasT;
            FrgCommon::frequency().interpolateOffset(s, lowerS, upperS, biasS);
            FrgCommon::frequency().interpolateOffset(t, lowerT, upperT, biasT);
 
            accessBuffer.frequencyWeights[0] = (1 - biasT) * (1 - biasS);
            accessBuffer.frequencyOffsets[0] = _generateAccessBufferOffset(lowerS, lowerT, exactU, accessBuffer.signFlag[0]);
            accessBuffer.frequencyWeights[1] = (1 - biasT) * biasS;
            accessBuffer.frequencyOffsets[1] = _generateAccessBufferOffset(upperS, lowerT, exactU, accessBuffer.signFlag[1]);
            accessBuffer.frequencyWeights[2] = biasT * (1 - biasS);
            accessBuffer.frequencyOffsets[2] = _generateAccessBufferOffset(lowerS, upperT, exactU, accessBuffer.signFlag[2]);
            accessBuffer.frequencyWeights[3] = biasT * biasS;
            accessBuffer.frequencyOffsets[3] = _generateAccessBufferOffset(upperS, upperT, exactU, accessBuffer.signFlag[3]);
        }
 
        return accessBuffer;
    }
 
    XYZVertexTwoParticleAccessBuffer<8> generateAccessBuffer(float s, float t, float u) const
    {
        XYZVertexTwoParticleAccessBuffer<8> accessBuffer;
 
        //map to positive frequency sector
        if (s < 0 && u < 0)
        {
            s = -s;
            u = -u;
        }
        else
        {
            if (s < 0)
            {
                s = -s;
                accessBuffer.siteExchange = true;
            }
            else if (u < 0)
            {
                u = -u;
                accessBuffer.siteExchange = true;
            }
        }
        if (t < 0)
        {
            t = -t;
        }
 
        int lowerS, upperS;
        float biasS;
        int lowerT, upperT;
        float biasT;
        int lowerU, upperU;
        float biasU;
 
        FrgCommon::frequency().interpolateOffset(s, lowerS, upperS, biasS);
        FrgCommon::frequency().interpolateOffset(t, lowerT, upperT, biasT);
        FrgCommon::frequency().interpolateOffset(u, lowerU, upperU, biasU);
 
        accessBuffer.frequencyWeights[0] = (1 - biasT) * (1 - biasS) * (1 - biasU);
        accessBuffer.frequencyOffsets[0] = _generateAccessBufferOffset(lowerS, lowerT, lowerU, accessBuffer.signFlag[0]);
        accessBuffer.frequencyWeights[1] = (1 - biasT) * biasS * (1 - biasU);
        accessBuffer.frequencyOffsets[1] = _generateAccessBufferOffset(upperS, lowerT, lowerU, accessBuffer.signFlag[1]);
        accessBuffer.frequencyWeights[2] = biasT * (1 - biasS) * (1 - biasU);
        accessBuffer.frequencyOffsets[2] = _generateAccessBufferOffset(lowerS, upperT, lowerU, accessBuffer.signFlag[2]);
        accessBuffer.frequencyWeights[3] = biasT * biasS * (1 - biasU);
        accessBuffer.frequencyOffsets[3] = _generateAccessBufferOffset(upperS, upperT, lowerU, accessBuffer.signFlag[3]);
        accessBuffer.frequencyWeights[4] = (1 - biasT) * (1 - biasS) * biasU;
        accessBuffer.frequencyOffsets[4] = _generateAccessBufferOffset(lowerS, lowerT, upperU, accessBuffer.signFlag[4]);
        accessBuffer.frequencyWeights[5] = (1 - biasT) * biasS * biasU;
        accessBuffer.frequencyOffsets[5] = _generateAccessBufferOffset(upperS, lowerT, upperU, accessBuffer.signFlag[5]);
        accessBuffer.frequencyWeights[6] = biasT * (1 - biasS) * biasU;
        accessBuffer.frequencyOffsets[6] = _generateAccessBufferOffset(lowerS, upperT, upperU, accessBuffer.signFlag[6]);
        accessBuffer.frequencyWeights[7] = biasT * biasS * biasU;
        accessBuffer.frequencyOffsets[7] = _generateAccessBufferOffset(upperS, upperT, upperU, accessBuffer.signFlag[7]);
 
        return accessBuffer;
    }
 
    float _directAccessMapFrequencyExchange(const int siteOffset, const int sOffset, const int tOffset, const int uOffset, const SpinComponent symmetry) const
    {
        ASSERT(siteOffset >= 0 && siteOffset < FrgCommon::lattice().size);
        ASSERT(sOffset >= 0 && sOffset < FrgCommon::frequency().size);
        ASSERT(tOffset >= 0 && tOffset < FrgCommon::frequency().size);
        ASSERT(uOffset >= 0 && uOffset < FrgCommon::frequency().size);
 
        if (sOffset >= uOffset) return _directAccess(siteOffset, sOffset, tOffset, uOffset, symmetry);
        else 
        {
            if (symmetry != SpinComponent::None) return _directAccess(siteOffset, uOffset, tOffset, sOffset, symmetry);
            else return -_directAccess(siteOffset, uOffset, tOffset, sOffset, symmetry);
        }
    }
 
    float _directAccess(const int siteOffset, const int sOffset, const int tOffset, const int uOffset, const SpinComponent symmetry) const
    {
        ASSERT(siteOffset >= 0 && siteOffset < FrgCommon::lattice().size);
        ASSERT(sOffset >= 0 && sOffset < FrgCommon::frequency().size);
        ASSERT(tOffset >= 0 && tOffset < FrgCommon::frequency().size);
        ASSERT(uOffset >= 0 && uOffset < FrgCommon::frequency().size);
        ASSERT(sOffset >= uOffset);
 
        if (symmetry == SpinComponent::X) return _dataXX[_memoryStepLatticeT * (sOffset * (sOffset + 1) / 2 + uOffset) + _memoryStepLattice * tOffset + siteOffset];
        if (symmetry == SpinComponent::Y) return _dataYY[_memoryStepLatticeT * (sOffset * (sOffset + 1) / 2 + uOffset) + _memoryStepLattice * tOffset + siteOffset];
        if (symmetry == SpinComponent::Z) return _dataZZ[_memoryStepLatticeT * (sOffset * (sOffset + 1) / 2 + uOffset) + _memoryStepLattice * tOffset + siteOffset];
        else return *(_dataDD + _memoryStepLatticeT * (sOffset * (sOffset + 1) / 2 + uOffset) + _memoryStepLattice * tOffset + siteOffset);
    }
 
    int _generateAccessBufferOffset(const int sOffset, const int tOffset, const int uOffset, int &signFlag) const
    {
        ASSERT(sOffset >= 0 && sOffset < FrgCommon::frequency().size);
        ASSERT(tOffset >= 0 && tOffset < FrgCommon::frequency().size);
        ASSERT(uOffset >= 0 && uOffset < FrgCommon::frequency().size);
 
        if (sOffset < uOffset)
        {
            signFlag = -1;
            return _memoryStepLatticeT * (uOffset * (uOffset + 1) / 2 + sOffset) + tOffset * _memoryStepLattice;
        }
        else
        {
            signFlag = 1;
            return _memoryStepLatticeT * (sOffset * (sOffset + 1) / 2 + uOffset) + tOffset * _memoryStepLattice;
        }
    }
 
    int size; 
    int sizeFrequency; 
 
    float *_dataXX; 
    float *_dataYY; 
    float *_dataZZ; 
    float *_dataDD; 
    int _memoryStepLatticeT; 
    int _memoryStepLattice; 
};