PropertyExtractor.h

Go to the documentation of this file.

/* Copyright 2016 Kristofer Björnson
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */

#ifndef COM_DAFER45_TBTK_PROPERTY_EXTRACTOR_PROPERTY_EXTRACTOR
#define COM_DAFER45_TBTK_PROPERTY_EXTRACTOR_PROPERTY_EXTRACTOR

#include "TBTK/HoppingAmplitudeSet.h"
#include "TBTK/Index.h"
#include "TBTK/Property/AbstractProperty.h"
#include "TBTK/Property/Density.h"
#include "TBTK/Property/DOS.h"
#include "TBTK/Property/LDOS.h"
#include "TBTK/Property/Magnetization.h"
#include "TBTK/Property/SpinPolarizedLDOS.h"
#include "TBTK/Solver/Solver.h"

#include <complex>
//#include <initializer_list>

namespace TBTK{
namespace PropertyExtractor{

class PropertyExtractor{
public:
    PropertyExtractor();

    virtual ~PropertyExtractor();

    virtual void setEnergyWindow(
        double lowerBound,
        double upperBound,
        int energyResolution
    );

    virtual void setEnergyWindow(
        int lowerFermionicMatsubaraEnergyIndex,
        int upperFermionicMatsubaraEnergyIndex,
        int lowerBosonicMatsubaraEnergyIndex,
        int upperBosonicMatsubaraEnergyIndex
    );

    virtual void setEnergyInfinitesimal(double energyInfinitesimal);

    virtual Property::Density calculateDensity(
        Index pattern,
        Index ranges
    );

    virtual Property::Density calculateDensity(
        std::vector<Index> patterns
    );

    virtual Property::Magnetization calculateMagnetization(
        Index pattern,
        Index ranges
    );

    virtual Property::Magnetization calculateMagnetization(
        std::vector<Index> patterns
    );

    virtual Property::LDOS calculateLDOS(Index pattern, Index ranges);

    virtual Property::LDOS calculateLDOS(
        std::vector<Index> patterns
    );

    virtual Property::SpinPolarizedLDOS calculateSpinPolarizedLDOS(
        Index pattern,
        Index ranges
    );

    virtual Property::SpinPolarizedLDOS calculateSpinPolarizedLDOS(
        std::vector<Index> patterns
    );

    virtual std::complex<double> calculateExpectationValue(
        Index to,
        Index from
    );

    virtual Property::DOS calculateDOS();

    virtual Property::DOS sampleDOS(
        unsigned int numSamples,
        const std::vector<Index> &patterns = {},
        unsigned int seed = time(nullptr)
    );

    virtual double calculateEntropy();
protected:
    enum class EnergyType{Real, Matsubara};

    EnergyType getEnergyType() const;

    int getEnergyResolution() const;

    double getLowerBound() const;

    double getUpperBound() const;

    int getLowerFermionicMatsubaraEnergyIndex() const;

    int getUpperFermionicMatsubaraEnergyIndex() const;

    int getLowerBosonicMatsubaraEnergyIndex() const;

    int getUpperBosonicMatsubaraEnergyIndex() const;

    class Information{
    public:
        Information();

        void setSpinIndex(int spinIndex);

        int getSpinIndex() const;
    private:
        int spinIndex;
    };

    /*** Get the nergy infinitesimal.
     *
     *  @return The energy infinitesimal. */
    double getEnergyInfinitesimal() const;

    template<typename DataType>
    void calculate(
        void (*callback)(
            PropertyExtractor *cb_this,
            Property::Property &property,
            const Index &index,
            int offset,
            Information &information
        ),
        Property::AbstractProperty<DataType> &property,
        Index pattern,
        const Index &ranges,
        int currentOffset,
        int offsetMultiplier,
        Information &information
    );

    template<typename DataType>
    void calculate(
        void (*callback)(
            PropertyExtractor *cb_this,
            Property::Property &property,
            const Index &index,
            int offset,
            Information &information
        ),
        const IndexTree &allIndices,
        const IndexTree &memoryLayout,
        Property::AbstractProperty<DataType> &abstractProperty,
        Information &information
    );

    void ensureCompliantRanges(const Index &pattern, Index &ranges);

    std::vector<int> getLoopRanges(
        const Index &pattern,
        const Index &ranges
    );

    IndexTree generateIndexTree(
        std::vector<Index> patterns,
        const HoppingAmplitudeSet &hoppingAmplitudeSet,
        bool keepSummationWildcards,
        bool keepSpinWildcards
    );
private:
    EnergyType energyType;

    static constexpr int ENERGY_RESOLUTION = 1000;

    int energyResolution;

    static constexpr double LOWER_BOUND = -1.;

    double lowerBound;

    static constexpr double UPPER_BOUND = 1.;

    double upperBound;

    static constexpr int LOWER_FERMIONIC_MATSUBARA_ENERGY_INDEX = -1;

    int lowerFermionicMatsubaraEnergyIndex;

    static constexpr int UPPER_FERMIONIC_MATSUBARA_ENERGY_INDEX = 1;

    int upperFermionicMatsubaraEnergyIndex;

    static constexpr int LOWER_BOSONIC_MATSUBARA_ENERGY_INDEX = 0;

    int lowerBosonicMatsubaraEnergyIndex;

    static constexpr int UPPER_BOSONIC_MATSUBARA_ENERGY_INDEX = 0;

    int upperBosonicMatsubaraEnergyIndex;

    static constexpr double ENERGY_INFINITESIMAL = 1e-3;

    double energyInfinitesimal;

    virtual const Solver::Solver& getSolver() const;
};

inline PropertyExtractor::EnergyType PropertyExtractor::getEnergyType() const{
    return energyType;
}

inline int PropertyExtractor::getEnergyResolution() const{
    TBTKAssert(
        energyType == EnergyType::Real,
        "PropertyExtractor::PropertyExtractor::getEnergyResolution()",
        "The energy resolution cannot be accessed when the energy type"
        << " is Matsubara.",
        "Use PropertyExtractor::PropertyExtractor::setEnergyWindow()"
        << " with three arguments if real energies are wanted for the"
        << " PropertyExtractor."
    );

    return energyResolution;
}

inline double PropertyExtractor::getLowerBound() const{
    TBTKAssert(
        energyType == EnergyType::Real,
        "PropertyExtractor::PropertyExtractor::getLowerBound()",
        "The lower bound cannot be accessed when the energy type is"
        << " Matsubara.",
        "Use PropertyExtractor::PropertyExtractor::setEnergyWindow()"
        << " with three arguments if real energies are wanted for the"
        << " PropertyExtractor."
    );

    return lowerBound;
}

inline double PropertyExtractor::getUpperBound() const{
    TBTKAssert(
        energyType == EnergyType::Real,
        "PropertyExtractor::PropertyExtractor::getUpperBound()",
        "The upper bound cannot be accessed when the energy type is"
        << " Matsubara.",
        "Use PropertyExtractor::PropertyExtractor::setEnergyWindow()"
        << " with three arguments if real energies are wanted for the"
        << " PropertyExtractor."
    );

    return upperBound;
}

inline int PropertyExtractor::getLowerFermionicMatsubaraEnergyIndex() const{
    TBTKAssert(
        energyType == EnergyType::Matsubara,
        "PropertyExtractor::PropertyExtractor::getLowerFermionicMatsubaraEnergyIndex()",
        "The lower Fermionic Matsubara energy index cannot be accessed"
        << " when the energy type is real.",
        "Use PropertyExtractor::PropertyExtractor::setEnergyWindow()"
        << " with four arguments if Matsubara energies are wanted for"
        << " the PropertyExtractor."
    );

    return lowerFermionicMatsubaraEnergyIndex;
}

inline int PropertyExtractor::getUpperFermionicMatsubaraEnergyIndex() const{
    TBTKAssert(
        energyType == EnergyType::Matsubara,
        "PropertyExtractor::PropertyExtractor::getUpperFermionicMatsubaraEnergyIndex()",
        "The upper Fermionic Matsubara energy index cannot be accessed"
        << " when the energy type is real.",
        "Use PropertyExtractor::PropertyExtractor::setEnergyWindow()"
        << " with four arguments if Matsubara energies are wanted for"
        << " the PropertyExtractor."
    );

    return upperFermionicMatsubaraEnergyIndex;
}

inline int PropertyExtractor::getLowerBosonicMatsubaraEnergyIndex() const{
    TBTKAssert(
        energyType == EnergyType::Matsubara,
        "PropertyExtractor::PropertyExtractor::getLowerBosonicMatsubaraEnergyIndex()",
        "The lower Bosonic Matsubara energy index cannot be accessed"
        << " when the energy type is real.",
        "Use PropertyExtractor::PropertyExtractor::setEnergyWindow()"
        << " with four arguments if Matsubara energies are wanted for"
        << " the PropertyExtractor."
    );

    return lowerBosonicMatsubaraEnergyIndex;
}

inline int PropertyExtractor::getUpperBosonicMatsubaraEnergyIndex() const{
    TBTKAssert(
        energyType == EnergyType::Matsubara,
        "PropertyExtractor::PropertyExtractor::getUpperBosonicMatsubaraEnergyIndex()",
        "The upper Bosonic Matsubara energy index cannot be accessed"
        << " when the energy type is real.",
        "Use PropertyExtractor::PropertyExtractor::setEnergyWindow()"
        << " with four arguments if Matsubara energies are wanted for"
        << " the PropertyExtractor."
    );

    return upperBosonicMatsubaraEnergyIndex;
}

inline double PropertyExtractor::getEnergyInfinitesimal() const{
    return energyInfinitesimal;
}

template<typename DataType>
void PropertyExtractor::calculate(
    void (*callback)(
        PropertyExtractor *cb_this,
        Property::Property &property,
        const Index &index,
        int offset,
        Information &information
    ),
    Property::AbstractProperty<DataType> &property,
    Index pattern,
    const Index &ranges,
    int currentOffset,
    int offsetMultiplier,
    Information &information
){
    //Find the next specifier index.
    int currentSubindex = pattern.getSize()-1;
    for(; currentSubindex >= 0; currentSubindex--){
        if(pattern.at(currentSubindex) < 0)
            break;
    }

    if(currentSubindex == -1){
        //No further specifier index found. Call the callback.
        callback(this, property, pattern, currentOffset, information);
    }
    else{
        //Ensure that the specifier is valid for the Ranges format.
        TBTKAssert(
            pattern.at(currentSubindex).isSummationIndex()
            || pattern.at(currentSubindex).isRangeIndex(),
            "PropertyExtractor::calculate()",
            "Invalid specifier found at subindex "
            << currentSubindex << ".",
            "Did you mean IDX_SUM_ALL, IDX_X, IDX_Y, or IDX_Z?"
        );

        //Multiply the memory offset for non summation indices.
        int nextOffsetMultiplier = offsetMultiplier;
        if(!pattern.at(currentSubindex).isSummationIndex())
            nextOffsetMultiplier *= ranges.at(currentSubindex);

        //Set flag indicating whether the current subindex is a
        //summation index.
        bool isSumIndex = false;
        if(pattern.at(currentSubindex).isSummationIndex())
            isSumIndex = true;

        //Recurively call the calculate function with the specifier at
        //the current subindex replaced by each subindex value in the
        //corresponding range.
        for(int n = 0; n < ranges.at(currentSubindex); n++){
            pattern.at(currentSubindex) = n;
            calculate(
                callback,
                property,
                pattern,
                ranges,
                currentOffset,
                nextOffsetMultiplier,
                information
            );
            if(!isSumIndex)
                currentOffset += offsetMultiplier;
        }
    }
}

template<typename DataType>
void PropertyExtractor::calculate(
    void (*callback)(
        PropertyExtractor *cb_this,
        Property::Property &property,
        const Index &index,
        int offset,
        Information &information
    ),
    const IndexTree &allIndices,
    const IndexTree &memoryLayout,
    Property::AbstractProperty<DataType> &abstractProperty,
    Information &information
){
    for(
        IndexTree::ConstIterator iterator = allIndices.cbegin();
        iterator != allIndices.end();
        ++iterator
    ){
        const Index &index = *iterator;
        std::vector<unsigned int> spinIndices
            = memoryLayout.getSubindicesMatching(
                IDX_SPIN,
                index,
                IndexTree::SearchMode::MatchWildcards
            );
        if(spinIndices.size() != 0){
            TBTKAssert(
                spinIndices.size() == 1,
                "PropertyExtractor::calculate()",
                "Several spin indeces found.",
                "Use IDX_SPIN at most once per pattern to"
                << " indicate spin index."
            );
            information.setSpinIndex(spinIndices[0]);
        }

        callback(
            this,
            abstractProperty,
            index,
            abstractProperty.getOffset(index),
            information
        );
    }
}

inline void PropertyExtractor::Information::setSpinIndex(int spinIndex){
    this->spinIndex = spinIndex;
}

inline int PropertyExtractor::Information::getSpinIndex() const{
    return spinIndex;
}

inline const Solver::Solver& PropertyExtractor::getSolver() const{
    TBTKExit(
        "PropertyExtractor::getSolver()",
        "Missing implementation of PropertyExtractor::getSolver() in"
        " deriving class.",
        ""
    );
}

};  //End of namespace PropertyExtractor
};  //End of namespace TBTK

#endif