graph.h

Go to the documentation of this file.

// This file is part of MARTY.
// 
// MARTY 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.
// 
// MARTY 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 MARTY. If not, see <https://www.gnu.org/licenses/>.

#ifndef GRAPH_H_INCLUDED
#define GRAPH_H_INCLUDED

#include "quantumField.h"
#include "fermionFlow.h"

namespace mty {
class FeynruleMomentum;
class FeynmanDiagram;
class ConjugationList;
class Model;
}

namespace mty::wick {

class Node: public std::enable_shared_from_this<Node>{

    public:

    explicit
    Node(const mty::QuantumField* t_field): field(t_field){};

    Node(const Node& node) = default;

    Node(Node&& node) = default;

    bool isFree() const {
        return (not partner.lock());
    };

    static std::vector<std::shared_ptr<Node>> copyGraph(
            std::vector<std::shared_ptr<Node>> const& init);

    public:

    const mty::QuantumField *field;

    std::weak_ptr<Node> partner;
};

// Forward declaration of ConnectedComponent class
class ConnectedComponent;

class Vertex: public std::vector<std::shared_ptr<Node>> {

    public:

    Vertex();

    explicit
    Vertex(const std::vector<const mty::QuantumField*>& fields,
           size_t                                       id,
           bool                                         t_external=false);

    Vertex(const Vertex& other);

    bool isExternal() const;

    bool isFree() const;

    bool hasFreeNode() const;

    bool hasFreeNode(const mty::QuantumField* field) const;

    bool containsExplicit(const mty::QuantumField* field) const;

    bool isDegenerateWith(const Vertex& other) const;

    const mty::QuantumField* getFirstContractibleField() const;

    std::shared_ptr<Node> getFirstContractibleNode() const; 

    std::shared_ptr<Node> getFirstContractibleNode(
            const mty::QuantumField* field) const;

    std::vector<std::shared_ptr<Node>> getContractibleNodes() const;

    std::vector<std::shared_ptr<Node>> getContractibleNodes(
            const mty::QuantumField* field) const;

    void setExternal(bool t_external);
    
    int getDegeneracy(const mty::QuantumField* t_field) const;

    int getTotalDegeneracyFactor() const;

    Vertex& operator=(const Vertex& other);

    bool operator==(const Vertex& other) const;

    bool operator!=(const Vertex& other) const;

    friend
    std::ostream& operator<<(std::ostream& fout,
                             const Vertex& v);

    protected:

    bool external;

    size_t id;
};

class ConnectedComponent {

    IMPLEMENTS_STD_VECTOR(Vertex, vertices)

    public:

    ConnectedComponent();

    explicit
    ConnectedComponent(const Vertex& firstVertex);

    ConnectedComponent(const ConnectedComponent& other) = default;

    std::vector<Vertex> const &getVertices() const;

    void addVertex(const Vertex& newVertex);

    bool hasValenceLeft() const;

    const mty::QuantumField* getFirstContractibleField() const;

    int getFirstContractibleVertex() const;

    std::vector<int> getContractibleVertices(
            const mty::QuantumField* field) const;

    std::shared_ptr<Node> getFirstContractibleNode() const;

    std::vector<std::shared_ptr<Node>> getNodes() const;

    int connect(int indexVertex);

    bool connect(Vertex& other);

    friend
    std::ostream& operator<<(std::ostream&             fout,
                             const ConnectedComponent& c);

    private:

    std::vector<Vertex> vertices;
};

class Graph {

    public:

    using Expr_type = std::pair<csl::Expr, ConjugationList>;

    Graph();

    explicit
    Graph(const std::vector<mty::QuantumField>& field,
          std::map<csl::Tensor, size_t>& vertexIds,
          bool                           ruleMode);

    explicit
    Graph(const csl::Expr& expr,
          std::map<csl::Tensor, size_t>& vertexIds,
          bool                           ruleMode);

    Graph(Graph const& other,
          std::vector<std::shared_ptr<Node>> const& init);

    ~Graph(){};

    Graph(const Graph& other) = default;

    Graph& operator=(const Graph& other) = default;
    // Graph(Graph const& other)
    //     :sortedNodes(other.sortedNodes),
    //     factor(other.factor),
    //     symmetrizeExternalLegs(other.symmetrizeExternalLegs),
    //     init(other.init),
    //     initialOrder(other.initialOrder),
    //     connectedCompo(other.connectedCompo),
    //     extVertex(other.extVertex),
    //     intVertex(other.intVertex),
    //     independentVertices(other.independentVertices),
    //     symmetry(other.symmetry),
    //     fixedExpression(other.fixedExpression),
    //     expression(other.expression)
    // {
    //     ++countCopy;
    // }

    // Graph& operator=(Graph const& other) {
    //     sortedNodes = other.sortedNodes;
    //     factor = other.factor;
    //     symmetrizeExternalLegs = other.symmetrizeExternalLegs;
    //     init = other.init;
    //     initialOrder = other.initialOrder;
    //     connectedCompo = other.connectedCompo;
    //     extVertex = other.extVertex;
    //     intVertex = other.intVertex;
    //     independentVertices = other.independentVertices;
    //     symmetry = other.symmetry;
    //     fixedExpression = other.fixedExpression;
    //     expression = other.expression;
    //     ++countCopy;
    //     return *this;
    // }

    public: 

    Graph copy() const;

    int getFactor() const;

    int getTotalDegeneracyFactor() const;

    std::vector<Vertex> const &getVertices() const;

    void addFactor(int t_factor);

    void setFactor(int t_factor);

    void setSymmetrizeExternalLegs(bool t_symmetrize);

    void multiply(int t_factor);

    void multiply(csl::Expr const &t_factor);

    bool isFullyConnected() const;

    bool isPhysical() const;

    bool isValid() const;

    int getNLoops() const;

    int getFieldDimension() const;

    std::vector<Expr_type> getExpression() const;

    std::vector<Expr_type> getExpression(
            mty::FeynruleMomentum& witness
            ) const;

    std::vector<const mty::QuantumField*> getFields() const;

    std::vector<const mty::QuantumField*> getConnectedFieldsInVertex(
            const mty::QuantumField* init) const;

    const mty::QuantumField* getConnectedField(
            const mty::QuantumField* init) const;

    std::vector<std::shared_ptr<Node>>& getNodes();

    std::vector<std::shared_ptr<Node>>const& getNodes() const;

    bool empty() const;

    std::vector<std::shared_ptr<Graph>> contractionStep() const;

    bool compare(const Graph& other,
                 std::map<csl::Tensor, csl::Tensor>& mapping,
                 bool         fieldBlind = false) const;

    bool compare(const Graph& other,
                 bool         fieldBlind = false) const;

    void contractHighMass(
            mty::QuantumFieldParent const *field
            );

    bool operator==(const Graph& other) const;

    friend std::ostream& operator<<(std::ostream& fout,
                                    const Graph&  g);

    static 
    Vertex const *getVertexOf(
            std::shared_ptr<Node> const &node,
            std::vector<Vertex> const &vertices
            );

    static
    std::vector<std::shared_ptr<Node>> nextNodes(
            std::shared_ptr<Node> const &node,
            std::vector<Vertex> const &vertices
            );

    static
    int countExternalLegs(
            std::vector<csl::Tensor>::iterator first,
            std::vector<csl::Tensor>::iterator last,
            std::vector<Vertex>          const &vertices
            );

    static
    int walk(
            std::vector<csl::Tensor>::iterator first,
            std::vector<csl::Tensor>::iterator last,
            std::shared_ptr<Node>        const &node,
            std::vector<Vertex>          const &vertices
            );

    private:
    
    void isolateIndependantVertices();

    void removeIntVertex(int indexVertex);

    void initConnectedComponent();

    Expr_type getPartialExpression(
            std::vector<std::shared_ptr<Node>>&   nodes,
            const std::vector<mty::QuantumField>& initialOrder,
            mty::FeynruleMomentum& witnessMapping,
            csl::Expr const&            globalFactor) const;

    csl::ObjectSymmetry<const mty::QuantumField*>
        getSymmetry(const Vertex&            vertex,
                    const mty::QuantumField* field) const;

    std::vector<std::vector<std::shared_ptr<Node>>> symmetrize() const;

    static std::vector<std::shared_ptr<Node>> copyNodes(
            const std::vector<std::shared_ptr<Node>>& toCopy);

    static void applySymmetry(
            std::vector<std::shared_ptr<Node>>&             nodes,
            const csl::ObjectPermutation<const mty::QuantumField*>& permutation);

    std::vector<int> getContractibleExtVertices(
            const mty::QuantumField* field) const;

    std::vector<int> getContractibleIntVertices(
            const mty::QuantumField* field) const;

    std::vector<int> getContractibleConnectedVertices(
            const mty::QuantumField* field) const;

    bool contractExternal(const mty::QuantumField* field,
                          int                 indexVertex);

    bool contractInternal(const mty::QuantumField* field,
                          int                 indexVertex);

    bool contractConnected(const mty::QuantumField* field,
                           int                      indexVertex);

    static void addFoundNode(const std::shared_ptr<Node>&      newNode,
                             std::vector<csl::Tensor>& foundNodes);

    static void sortNodes(std::vector<std::shared_ptr<Node>>& nodes);

    static bool compareFieldsDummyPoint(
            const mty::QuantumField*                            fieldA,
            const mty::QuantumField*                            fieldB,
            std::map<csl::Tensor, csl::Tensor>& constraints,
            bool fieldBlind = false);

    static bool compareNodesWithConstraints(
            const std::shared_ptr<Node>&                        nodeA,
            const std::shared_ptr<Node>&                        nodeB,
            std::map<csl::Tensor, csl::Tensor>& constraints,
            bool fieldBlind = false);

    private:

    mutable
    std::vector<std::shared_ptr<Node>> sortedNodes;
    long long int factor = 1;

    csl::Expr          exprFactor = CSL_1;

    bool symmetrizeExternalLegs = false;

    bool ruleMode;

    std::shared_ptr<std::vector<Vertex>> init;

    std::shared_ptr<std::vector<mty::QuantumField>> initialOrder;

    ConnectedComponent  connectedCompo;

    std::vector<Vertex> extVertex;

    std::vector<Vertex> intVertex;

    std::vector<int>    independentVertices;

    csl::ObjectSymmetry<const mty::QuantumField*> symmetry 
        = Id_Sym(const mty::QuantumField*);

    mutable bool fixedExpression = false;

    mutable std::vector<Expr_type> expression;
};

class WickCalculator {

    public:

    static csl::Expr applyWickTheoremOnDiagram(
            const Graph& diagram,
            std::vector<mty::FeynruleMomentum>& witnessMapping,
            bool ruleMode = true);

    static csl::vector_expr applyWickTheoremOnDiagrams(
            const std::vector<std::shared_ptr<Graph>>& diagrams,
            std::vector<mty::FeynruleMomentum>& witnessMapping,
            bool ruleMode = true);

    static std::vector<mty::FeynmanDiagram> getDiagrams(
            mty::Model const *model,
            const csl::Expr& initialExpr,
            std::map<csl::Tensor, size_t>& vertexIds,
            std::vector<mty::FeynruleMomentum>& witnessMapping,
            bool symmetriseExternalLegs = false,
            bool ruleMode = true);

    static std::vector<mty::FeynmanDiagram> getDiagrams(
            mty::Model const *model,
            const csl::Expr& initialExpr,
            std::map<csl::Tensor, size_t>& vertexIds,
            bool symmetriseExternalLegs = false,
            bool ruleMode = true);

    static bool isContractionZero(
            const std::vector<mty::QuantumField>& fields,
            size_t maxLoops  = size_t(-1),
            size_t nVertices = size_t(-1));

    static bool isContractionZero(
            const std::vector<mty::QuantumField const*>& fields,
            size_t maxLoops  = size_t(-1),
            size_t nVertices = size_t(-1));

    private:

    static std::vector<std::vector<mty::QuantumField>> splitFields(
            const std::vector<mty::QuantumField>& field);

    static std::vector<std::vector<mty::QuantumField const*>> splitFields(
            const std::vector<mty::QuantumField const*>& field);

    private:

    WickCalculator();

    explicit
    WickCalculator(const Graph& freeDiagram);

    WickCalculator(const WickCalculator& other) = default;

    ~WickCalculator(){};

    void setSymmetrizeExternalLegs(bool t_symmetrize);

    Graph getInitialDiagram() const;

    std::vector<std::shared_ptr<Graph>> getDiagrams();

    void eliminateNonPhysicalDiagrams(
            std::vector<std::shared_ptr<Graph>>& diagrams);

    void calculateDiagrams();

    private:

    bool diagramsCalculated = false;

    bool symmetrizeExternalLegs = false;

    Graph initialDiagram;

    std::vector<std::shared_ptr<Graph>> feynmanDiagram;
};

/*************************************************/
// Functions                                     //
/*************************************************/

void contractNodes(const std::shared_ptr<Node>& node1,
                   const std::shared_ptr<Node>& node2);

bool areDegenerate(const mty::QuantumField& f1, 
                   const mty::QuantumField& f2);

bool areExactlyContractible(const mty::QuantumField& f1, 
                            const mty::QuantumField& f2);

std::ostream& operator<<(std::ostream&                fout,
                         const std::shared_ptr<Node>& node);

bool operator==(const std::shared_ptr<Node>& A,
                const std::shared_ptr<Node>& B);

bool comparePriority(const std::shared_ptr<Node>&            A,
                     const std::shared_ptr<Node>&            B,
                     const std::vector<csl::Tensor>& foundNodes);

bool internal_comparePriority(const std::shared_ptr<Node>&            A,
                              const std::shared_ptr<Node>&            B,
                              const std::vector<csl::Tensor>& foundNodes);

std::vector<mty::QuantumField> convertExprToFields(const csl::Expr& expr);

int countFermions(std::vector<const mty::QuantumField*>::const_iterator begin,
                  std::vector<const mty::QuantumField*>::const_iterator end);

int getCommutationSign(const std::vector<mty::QuantumField>& A,
                       std::vector<const mty::QuantumField*>       B);

int getCommutationSign(const std::vector<const mty::QuantumField*>& A,
                       std::vector<const mty::QuantumField*>       B);

csl::vector_expr convertGraphsToCorrelators(
        const std::vector<std::shared_ptr<Graph>>& diagrams,
        std::vector<mty::FeynruleMomentum>& witnessMapping,
        bool ruleMode = true); 

} // End of namespace mty::wick

#endif