#ifndef _ULLMANN_H
#define _ULLMANN_H

#include "Graph.h"

class Ullmann
{
    Graph &g, &q;
    VertexType qn, gn;
    std::vector<VertexType> q2g;                    // M
    std::vector<std::vector<VertexType>> candi_q2g; // M'

    const std::string &output_fname;
    int mapping_num;

public:
    Ullmann(Graph &g, Graph &q, const std::string &output_fname) : g(g), q(q), qn(q.N), gn(g.N), output_fname(output_fname), mapping_num(0) {}

    bool NLF()
    {
        VertexType qn = q.N, gn = g.N;
        candi_q2g.resize(qn);
        for (VertexType u = 0; u < qn; ++u)
        {
            VertexLabelType ulb = q.vertex_label[u];
            for (VertexType v = 0; v < gn; ++v)
            {
                if (g.vertex_label[v] == ulb && q.nbr_label_freq_cover(u, g.nbr_label_freq[v]))
                    candi_q2g[u].push_back(v);
            }
            if (candi_q2g[u].empty())
                return false;
        }
        return true;
    }
    bool filter(VertexType u, VertexType v)
    {
        // matchings of all nbrs of u before u, must be nbrs of v
        for (VertexType un = 0; un < u; ++un)
        {
            VertexType vn = q2g[un];
            if (q.mat[u * qn + un] != g.mat[v * gn + vn])
                return false;
        }
        return true;
    }
    bool check_q2g()
    {
        // M*MB: each entry in q2g: q2g[u]'s nbrs in g
        // check each 1 in MA with M*(M*MB)T
        size_t pos = 0;
        for (VertexType u = 0; u < qn; ++u)
        {
            for (VertexType unbr = 0; unbr < qn; ++unbr, ++pos)
            {
                EdgeLabelType uelb = q.mat[pos];
                if (uelb >= 0)
                {
                    // check MA[u][unbr]
                    VertexType v = q2g[unbr];           // (M*MB)T[unbr] is MB[v]
                    if (g.mat[v * gn + q2g[u]] != uelb) // M only has only one non-zero, thus check MB[v][q2g[u]] is enough
                        return false;
                }
            }
            // hey you see that
            // this is just checking whether each edge in q is map to an edge in g
        }
        return true;
    }
    void output_q2g()
    {
#ifndef SILENCE_RESULT_OUTPUT
        std::ofstream fout(output_fname, std::ios::app);
        fout << "# " << mapping_num++ << std::endl;
        for (VertexType u = 0; u < qn; ++u)
            fout << u << " " << q2g[u] << std::endl;
#endif
    }
    void init_dfs()
    {
        q2g.resize(q.N);
    }
    void dfs(VertexType u)
    {
        // std::cout << u << std::endl;
        if (u >= qn)
        {
            if (check_q2g())
                output_q2g();
            return;
        }
        for (VertexType v : candi_q2g[u])
        {
            if (filter(u, v) == false)
                continue;
            q2g[u] = v;
            dfs(u + 1);
        }
    }

    void run()
    {
        NLF();
        init_dfs();
        dfs(0);
    }
};

#endif // _ULLMANN_H