semiSimpleAlgebra.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 SEMISIMPLEALGEBRA_H_INCLUDED
#define SEMISIMPLEALGEBRA_H_INCLUDED

#include <vector>
#include <algorithm>
#include <iostream>
#include <initializer_list>
#include <string>
#include <csl.h>
#include "algebraState.h"

//Forward declaration of Irrep and SumIrrep
namespace mty {
class Irrep;
class SumIrrep;
}

namespace mty::algebra {

    enum class Type {
        R, 
        A, 
        B, 
        C, 
        D, 
        E6, 
        E7, 
        E8, 
        F4, 
        G2, 
    };

    std::ostream& operator<<(std::ostream& fout, mty::algebra::Type obj);
}

namespace mty {

int getCoxeterNumber(
        algebra::Type type,
        int           l = 0
        );

struct IndexData {
    AlgebraState state;
    size_t p;
    csl::Expr index;
};

class SemiSimpleAlgebra{

    public:

    explicit
    SemiSimpleAlgebra(int t_l);

    virtual
    ~SemiSimpleAlgebra(){};

    csl::Expr getQuadraticCasimir(const AlgebraState& irrep) const;

    virtual csl::Expr getIndex(
            const Irrep& irrep,
            size_t       n) const;

    virtual
    Irrep highestWeightRep(const AlgebraState& highestWeight, bool mult = true) const;

    virtual
    Irrep highestWeightRep(const std::vector<int>& highestWeight, bool mult = true) const;

    virtual
    Irrep getTrivialRep() const;

    virtual 
    Irrep getDefiningRep() const;

    virtual
    SumIrrep tensorProduct(const Irrep& A,
                           const Irrep& B,
                           bool mult = true) const; 
    virtual
    SumIrrep tensorProduct(const SumIrrep& A,
                           const Irrep& B,
                           bool mult = true) const; 
    virtual
    SumIrrep tensorProduct(const Irrep& A,
                           const SumIrrep& B,
                           bool mult = true) const; 
    SumIrrep tensorProduct(const SumIrrep& A,
                           const SumIrrep& B,
                           bool mult = true) const; 

    virtual
    void sortRep(std::vector<AlgebraState>& rep,
                 std::vector<int>&          mult) const;

    int getOrderL() const;

    virtual
    algebra::Type getType() const = 0;

    protected:

    void init();

    virtual
    void setCartanMatrix() = 0;

    virtual
    void setSquaredNorm() = 0;

    virtual csl::Expr getEta() const {
        return CSL_2;
    }

    private:

    static bool isPositiveRoot(const std::vector<int>& root);

    static bool isPositiveWeight(const AlgebraState& state);

    static AlgebraState addWeights(const AlgebraState& A,
                                   const AlgebraState& B);

    bool isOnWall(const AlgebraState& state) const;

    AlgebraState toDominantWeylChamber(const AlgebraState& state,
                                       int&                sign) const;

    std::vector<int> weylReflection(int                     simpleRoot,
                                    const std::vector<int>& root) const;

    AlgebraState weylReflection(int                 simpleRoot,
                                const AlgebraState& state) const;

    std::vector<int> weylReflection(int simpleRoot1,
                                    int simpleRoot2) const;

    void invertCartanMatrix();

    void generateWeylGroup();

    csl::Expr convertState(const AlgebraState& state) const;

    bool comparatorHighest(const AlgebraState& A,
                           const AlgebraState& B) const;

    csl::Expr computeScalarDot(const csl::Expr& rootA,
                          const csl::Expr& rootB) const;

    csl::Expr computeSquaredNorm(const csl::Expr& root) const;

    AlgebraState applyAnnihilationOperator(int                 annihilator,
                                           const AlgebraState& state) const;

    AlgebraState applyAnnihilationOperator(const std::vector<int>& root,
                                           const AlgebraState&     state) const;

    AlgebraState applyCreationOperator(int                 creator,
                                       const AlgebraState& state) const;

    AlgebraState applyCreationOperator(const std::vector<int>& root,
                                       const AlgebraState&     state) const;

    std::vector<AlgebraState> getSingleChain(
            const AlgebraState& state,
            int                 direction) const;

    void getSingleChainExperimental(
            const AlgebraState& state,
            std::vector<AlgebraState> &states,
            std::vector<std::vector<bool>> &directionExplored,
            int                 direction) const;

    std::vector<AlgebraState> getRootChain(
            const AlgebraState&             highestWeightState,
            std::vector<AlgebraState>       states 
                = std::vector<AlgebraState>(0),
            std::vector<std::vector<bool> > directions
                = std::vector<std::vector<bool> >(0)) const;

    void getRootChainExperimental(
            const AlgebraState&              highestWeightState,
            std::vector<AlgebraState>       &states,
            std::vector<std::vector<bool> > &directions) const;

    void computeMultiplicity(std::vector<int>&                multiplicities,   
                             const std::vector<AlgebraState>& rep,
                             const AlgebraState&              state) const;

    std::vector<int> multiplicities(std::vector<AlgebraState>& rep) const;

    friend
    std::ostream& operator<<(std::ostream&            fout,
                             const SemiSimpleAlgebra& algebra);

    protected:

    const int l;

    std::vector<int> rho;

    csl::vector_expr absoluteNorms;

    std::vector<int> squaredNorm;

    std::vector<std::vector<int>> positiveRoots;

    std::vector<std::vector<int>> cartanMatrix;

    std::vector<std::vector<int>> rootMetric;

    csl::Expr symbolicRho;

    csl::Expr symbolicCartan;

    csl::Expr inverseCartan;

    csl::Expr symbolicMetric;

    csl::vector_expr symbolicRoots;

    csl::vector_expr symbolicSquaredNorms;

    mutable
    std::vector<std::pair<AlgebraState, Irrep>> irreps;

    mutable std::vector<IndexData> indices;
};

} // End of namespace mty

namespace mty::algebra{

class R: public mty::SemiSimpleAlgebra{

    public:

    R();

    ~R(){};

    algebra::Type getType() const override;

    void setCartanMatrix() override;

    void setSquaredNorm() override;

    mty::Irrep highestWeightRep(
            const mty::AlgebraState& highestWeight,
            bool mult = true) const override;

    mty::Irrep highestWeightRep(
            const std::vector<int>& highestWeight,
            bool mult = true) const override;

    mty::Irrep getTrivialRep() const override;

    mty::SumIrrep tensorProduct(const mty::Irrep& A,
                                const mty::Irrep& B,
                                bool mult = true) const override;

    void sortRep(std::vector<mty::AlgebraState>& rep,
                 std::vector<int>&                mult) const override;

    csl::Expr getEta() const override {
        return CSL_0;
    }
};

class A: public mty::SemiSimpleAlgebra{

    public:

    explicit
    A(int t_l);

    ~A(){};

    algebra::Type getType() const override;

    void setCartanMatrix() override;

    void setSquaredNorm() override;

    csl::Expr getEta() const override {
        return CSL_1;
    }

    mty::Irrep getDefiningRep() const override;
};

class B: public mty::SemiSimpleAlgebra{

    public:

    explicit
    B(int t_l);

    ~B(){};

    algebra::Type getType() const override;

    void setCartanMatrix() override;

    void setSquaredNorm() override;

    mty::Irrep getDefiningRep() const override;
};

class C: public mty::SemiSimpleAlgebra{

    public:

    explicit
    C(int t_l);

    ~C(){};

    algebra::Type getType() const override;

    void setCartanMatrix() override;

    void setSquaredNorm() override;

    mty::Irrep getDefiningRep() const override;
};

class D: public mty::SemiSimpleAlgebra{

    public:

    explicit
    D(int t_l);

    ~D(){};

    algebra::Type getType() const override;

    void setCartanMatrix() override;

    void setSquaredNorm() override;

    mty::Irrep getDefiningRep() const override;
};

class E6: public mty::SemiSimpleAlgebra{

    public:

    E6();

    ~E6(){};

    algebra::Type getType() const override;

    void setCartanMatrix() override;

    void setSquaredNorm() override;

    mty::Irrep getDefiningRep() const override;
};

class E7: public mty::SemiSimpleAlgebra{

    public:

    E7();

    ~E7(){};

    algebra::Type getType() const override;

    void setCartanMatrix() override;

    void setSquaredNorm() override;

    mty::Irrep getDefiningRep() const override;
};

class E8: public mty::SemiSimpleAlgebra{

    public:

    E8();

    ~E8(){};

    algebra::Type getType() const override;

    void setCartanMatrix() override;

    void setSquaredNorm() override;

    mty::Irrep getDefiningRep() const override;
};

class F4: public mty::SemiSimpleAlgebra{

    public:

    F4();

    ~F4(){};

    algebra::Type getType() const override;

    void setCartanMatrix() override;

    void setSquaredNorm() override;

    mty::Irrep getDefiningRep() const override;
};

class G2: public mty::SemiSimpleAlgebra{

    public:

    G2();

    ~G2(){};

    algebra::Type getType() const override;

    void setCartanMatrix() override;

    void setSquaredNorm() override;

    mty::Irrep getDefiningRep() const override;
};
} // End of namespace mty::algebra

/*************************************************/
// Helper functions in calculations in .cpp      //
/*************************************************/

template<typename T>
std::ostream& operator<<(std::ostream& fout, const std::vector<T>& vec)
{
    for (const auto& el : vec)
        fout << el << "  ";
    fout << std::endl;

    return fout;
}

template<typename T>
int findInVector(const std::vector<T>& v,
                 const T             & element)
{
    const typename std::vector<T>::const_iterator iter
        = find(v.begin(), v.end(), element);
    if (iter != v.end())
        return distance(v.begin(), iter);
    return -1;
}

template<typename T>
int binaryFindInVector(const std::vector<T>& v,
                       const T             & element)
{
    size_t a = 0;
    size_t b = v.size();
    while (b != a) {
        const auto middleIndex = (a+b) / 2;
        const auto middle = v[middleIndex];
        if (element < middle) {
            b = middleIndex;
        }
        else if (middle < element) {
            a = middleIndex;
        }
        else 
            return middleIndex;
    }
    return -1;
}

template<typename T>
bool addCarefully(std::vector<T>& vec,
                 T                element)
{
    if (findInVector<T>(vec,element) == -1) {
        vec.push_back(element);
        return true;
    }
    return false;
}

template<typename T>
std::vector<T> operator+(const std::vector<T>& x,
                         const std::vector<T>& y)
{
    if (x.size() != y.size()) {
        std::cout << "Type: two std::vectors must have "
                  << "the same dimension in product.\n";
        std::exit(1);
    }
    std::vector<T> rep(x.size());
    for (size_t i=0; i!=x.size(); ++i)
        rep[i] = x[i]+y[i];

    return rep;
}

template<typename T>
std::vector<T> operator-(const std::vector<T>& x,
                         const std::vector<T>& y)
{
    if (x.size() != y.size()) {
        std::cout << "Type: two std::vectors must have "
                  << "the same dimension in product.\n";
        std::exit(1);
    }
    std::vector<T> rep(x.size());
    for (size_t i=0; i!=x.size(); ++i)
        rep[i] = x[i]-y[i];

    return rep;
}

template<typename T>
std::vector<T> operator-(const std::vector<T>& x)
{
    std::vector<T> rep(x);
    for (auto& el : rep)
        el *= -1;
    return rep;
}

template<typename T>
std::vector<T> operator*(T                     a,
                         const std::vector<T>& x)
{
    std::vector<T> rep(x);
    for (auto& el : rep)
        el *= a;
    return rep;
}

template<typename T>
std::vector<T> operator*(const std::vector<T>& x,
                         T                     a)
{
    std::vector<T> rep(x);
    for (auto& el : rep)
        el *= a;
    return rep;
}

template<typename T>
T operator*(const std::vector<T>& x,
            const std::vector<T>& y)
{
    if (x.size() != y.size()) {
        std::cout << "Type: two std::vectors must have "
                  << "the same dimension in product.\n";
        std::exit(1);
    }
    T rep = 0;
    for (size_t i=0; i!=x.size(); ++i)
        rep += x[i]*y[i];

    return rep;
}


template<typename T>
std::vector<T> operator*(const std::vector<std::vector<T>>& A,
                         const std::vector<T>&              x)
{
    if (A.size() == 0)
        return std::vector<T>(0);
    if (A[0].size() != x.size()) {
        std::cout << "Type: two std::vectors must have "
                  << "the same dimension in product.\n";
        std::exit(1);
    }
    std::vector<T> rep(A.size(),0);
    for (size_t i=0; i!=A.size(); ++i)
        for (size_t j=0; j!=x.size(); ++j)
            rep[i] += A[i][j]*x[j];

    return rep;
}

template<typename T>
std::vector<T> operator*(const std::vector<T>              &x,
                         const std::vector<std::vector<T>>& A)
{
    if (A.size() == 0)
        return std::vector<T>(0);
    if (A[0].size() != x.size()) {
        std::cout << "Type: two std::vectors must have "
                  << "the same dimension in product.\n";
        std::exit(1);
    }
    std::vector<T> rep(A.size(),0);
    for (size_t i=0; i!=A.size(); ++i)
        for (size_t j=0; j!=x.size(); ++j)
            rep[i] += A[i][j]*x[j];

    return rep;
}

template<typename T>
std::vector<std::vector<T> > operator*(const std::vector<std::vector<T>>& A,
                                       const std::vector<std::vector<T>>& B)
{
    if (A.size() == 0 or B.size() == 0) {
        return std::vector<std::vector<T> >(0);
    }
    if (A[0].size() != B.size()) {
        std::cout << "Type: two std::vectors must have "
                  << "the same dimension in product.\n";
        std::exit(1);
    }
    std::vector<std::vector<T> > rep(A.size(), std::vector<T>(B[0].size(),0));
    for (size_t i=0; i!=A.size(); ++i)
        for (size_t j=0; j!=B[0].size(); ++j)
            for (size_t k=0; k!=A[0].size(); ++k)
                rep[i][j] += A[i][k]*B[k][j];

    return rep;
}

#endif