https://github.com/uhebertj/chordal
Raw File
Tip revision: c0c3883cf637a6dd6b5a601f059881a01e39f4f5 authored by Ursula Hebert-Johnson on 15 July 2022, 07:45:27 UTC
treewidth -> max clique size
Tip revision: c0c3883
chordal.cpp
/*
 * DP algorithm to count labeled chordal graphs.
 */

#include<iostream>
#include<vector>
#include<string>
#include <boost/multiprecision/cpp_int.hpp>

using namespace boost::multiprecision;
using namespace std;

typedef int1024_t lli;

bool verbatim = false;
int rec_depth = 0;
lli notSaved = -1;
int w;

lli choose(int n, int k);

lli chordal(int k);
lli comp_chordal(int k);
vector<lli> table_chordal;

lli chordal_conn(int k);
lli comp_chordal_conn(int k);
vector<lli> table_chordal_conn;

lli g(int t, int x, int z, int k);
lli comp_g(int t, int x, int z, int k);
vector<vector<vector<vector<lli> > > > table_g;

lli gTilde(int t, int x, int z, int k);
lli comp_gTilde(int t, int x, int z, int k);
vector<vector<vector<vector <lli> > > > table_gTilde;

lli gHat(int t, int x, int z, int k);
lli comp_gHat(int t, int x, int z, int k);
vector<vector<vector<vector <lli> > > > table_gHat;

lli g1Tilde(int t, int x, int k);
lli comp_g1Tilde(int t, int x, int k);
vector<vector<vector<lli> > > table_g1Tilde;

lli g2Tilde(int t, int x, int k);
lli comp_g2Tilde(int t, int x, int k);
vector<vector<vector<lli> > > table_g2Tilde;

lli f(int t, int x, int l, int k);
lli comp_f(int t, int x, int l, int k);
vector<vector<vector<vector <lli> > > > table_f;

lli fTilde(int t, int x, int l, int k);
lli comp_fTilde(int t, int x, int l, int k);
vector<vector<vector<vector <lli> > > > table_fTilde;

lli fHat(int t, int x, int l, int k);
lli comp_fHat(int t, int x, int l, int k);
vector<vector<vector<vector <lli> > > > table_fHat;

lli fHat(int t, int x, int z, int l, int k);
lli comp_fHat(int t, int x, int z, int l, int k);
vector<vector<vector<vector <vector<lli> > > > > table_fHat5;

vector<vector<lli> > table_choose;


// DEBUG PRINTOUTS
void debugIn(string fname, lli x1, lli x2, lli x3) {
	cout << string(rec_depth, ' ') << fname << "(" << x1 << "," << x2 <<	"," << x3 << ")" << endl;}

void debugIn(string fname, lli x1, lli x2, lli x3, lli x4) {
	cout << string(rec_depth, ' ') << fname << "(" << x1 << "," << x2 <<	"," << x3 << "," << x4 << ")" << endl;}

void debugIn(string fname, lli x1, lli x2, lli x3, lli x4, lli x5) {
	cout << string(rec_depth, ' ') << fname << "(" << x1 << "," << x2 <<	"," << x3 << "," << x4 << "," << x5 <<")" << endl;}

void debugOut(string fname, lli x1, lli x2, lli x3, lli ans) {
	cout << string(rec_depth, ' ') << fname << "(" << x1 << "," << x2 <<	"," << x3 << ") = " << ans << endl;}

void debugOut(string fname, lli x1, lli x2, lli x3, lli x4, lli ans) {
	cout << string(rec_depth, ' ') << fname << "(" << x1 << "," << x2 <<	"," << x3 << "," << x4 << ") = " << ans << endl;}

void debugOut(string fname, lli x1, lli x2, lli x3, lli x4, lli x5, lli ans) {
	cout << string(rec_depth, ' ') << fname << "(" << x1 << "," << x2 <<	"," << x3 << "," << x4 << "," << x5 << ") = " << ans << endl;}



// fill table for binomial coefficients.
void chooseInit(int n) {
	table_choose.resize(n+1, vector<lli> (n+1, 0));
	table_choose[0][0] = 1;
	
	for(int m = 1; m <= n; ++m) {
		table_choose[m][0] = 1;
		for(int k = 1; k <= m; ++k) {
			table_choose[m][k] = table_choose[m-1][k] + table_choose[m-1][k-1];
		}
	}
}

// 10 choose 3 = 10*9*8 / 1 * 2 * 3
// n choose k for n > k defaults to 0.
lli choose(int n, int k) {
	if(n < k) return 0;
	return table_choose[n][k];
}

// Counts connected chordal graphs on k vertices
// assume k > 0
lli chordal_conn(int k) {
	rec_depth++;
	
	if (table_chordal_conn[k] == notSaved) table_chordal_conn[k] = comp_chordal_conn(k);
	lli ans = table_chordal_conn[k];
	
	rec_depth--;

	return ans;
}

lli comp_chordal_conn(int k) {
	//empty graph counts here
	// clique evap at time 1 countrs here for t>=2, not for t=1
	lli ans = 0;
		
	// clique evap at time 1 countrs here for t=1, not higher t. 
	for(int t = 1; t <= k; ++t) {
		for(int l = 1; l <= k; ++l) {
			ans += choose(k,l) * f(t, 0, l, k-l);
		}
	}

	return ans;
}

// Counts chordal graphs on k vertices
// assume k > 0
lli chordal(int k) {		
	rec_depth++;
	
	if (table_chordal[k] == notSaved) table_chordal[k] = comp_chordal(k);
	lli ans = table_chordal[k];
	
	rec_depth--;

	return ans;
}

lli comp_chordal(int k) {
	// Base case k = 0
	if(k == 0) return 1;

	// Guess the set of vertices in the first component
	lli ans = 0;
	for (int kk = 1; kk <= k; ++kk) {
		ans += choose(k - 1, kk - 1) * chordal_conn(kk) * chordal(k - kk);
	}
	
	return ans;
}

// chordal graphs that evaporate at or before t
// x special vertices (do NOT allow 0)
// k normal vertices
lli g(int t, int x, int z, int k) {
	
	if(verbatim) debugIn("g",t,x,z,k);
		
	rec_depth++;
	
	if(table_g[t][x][z][k] == notSaved) table_g[t][x][z][k] = comp_g(t, x, z, k);
	lli ans = table_g[t][x][z][k];
	
	rec_depth--;
	
	if(verbatim) debugOut("g",t,x,z,k,ans);
	
	return ans;

}

lli comp_g(int t, int x, int z, int k) {
	// Base case t=0
	if(t == 0) return k == 0 ? 1 : 0;
		
	// non zero t. If x is zero cant use x for connectivity, so have
	// to ensure connectivity some other way. Call f maybe??
	if(x == 0) {
		// EXCEPTION HERE
		exit(1);
	}	
		
	// main case: t > 0 and x > 0
	// Guess the set of vertices in components evaporate at time
	// Precisely t.	
	lli ans = 0;
	for(int kk = 0; kk <= k; ++kk) {
		ans += choose(k, kk) * gTilde(t, x, z, kk) * g(t - 1, x, z, k - kk);
	}
	
	return ans;
}

// Same as g, but now all connected components evaporate at time exactly t
// graphs where V = X count (vacously true)
//
// Disallowing x=0
//
// Guess component that contains lowest labeled vertex.
lli gTilde(int t, int x, int z, int k) {
	
	if(verbatim) debugIn("gTilde", t, x, z, k);
	
	rec_depth++;
	
	if(table_gTilde[t][x][z][k] == notSaved) table_gTilde[t][x][z][k] = comp_gTilde(t, x, z, k);
	lli ans = table_gTilde[t][x][z][k];
	
	rec_depth--;
	
	if(verbatim) debugOut("gTilde",t,x,z,k,ans);
	
	return ans;
	
}

lli comp_gTilde(int t, int x, int z, int k) {
	// base case: k=0
	if(k == 0) return 1;
	
	// if x=0 we cant use X for connectivity. Either define this
	// or remove x=0 from domain 
	if(x == 0) {
		// EXCEPTION HERE
		exit(1);
	}
	
	lli ans = 0;
	
	// t and x both nonzero.
	// guess #labels in comp.
	// guess its neigh into x
	for(int kk = 1; kk <= k; ++kk) {
		for(int xx = 1; xx <= x; ++xx) {
			ans += choose(k-1, kk-1) * (choose(x, xx) - choose(z, xx)) * g1Tilde(t, xx, kk) * gTilde(t, x, z, k-kk);
		}
	}
	
	return ans;
	
}

// Same as gTilde, but now no component of G-X sees all of X
// graphs where V = X count (vacously true)
//
// Disallowing x=0
//
// Guess component that contains lowest labeled vertex.
lli gHat(int t, int x, int z, int k) {
	
	if(verbatim) debugIn("gHat",t,x,z,k);
	
	rec_depth++;
	
	if(table_gHat[t][x][z][k] == notSaved) table_gHat[t][x][z][k] = comp_gHat(t, x, z, k);
	lli ans = table_gHat[t][x][z][k];

	rec_depth--;
	
	if(verbatim) debugOut("gHat",t,x,z,k,ans);
	
	return ans;
}

lli comp_gHat(int t, int x, int z, int k) {

	// base case: k=0
	if(k == 0) return 1;
	
	// if x=0 we cant use X for connectivity. Either define this
	// or remove x=0 from domain 
	if(x == 0) {
		// EXCEPTION HERE
		exit(1);
	}
	
	lli ans = 0;
	
	// t and x both nonzero.
	// guess #labels in comp.
	// guess its neigh into x
	// now xx < x because no component sees all of x.
	for(int kk = 1; kk <= k; ++kk) {
		for(int xx = 1; xx < x; ++xx) {
			ans += choose(k-1, kk-1) * (choose(x, xx) - choose(z, xx)) * g1Tilde(t, xx, kk) * gHat(t, x, z, k-kk);
		}
	}
	
	return ans;
}

// exactly one comp of G-X, evap at exactly t. Sees all of X.
// Disallowing x=0 
lli g1Tilde(int t, int x, int k) {
	
	if(verbatim) debugIn("g1Tilde", t, x, k);
	
	rec_depth++;
	
	if(table_g1Tilde[t][x][k] == notSaved) table_g1Tilde[t][x][k] = comp_g1Tilde(t, x, k);
	lli ans = table_g1Tilde[t][x][k];
	
	rec_depth--;
	
	if(verbatim) debugOut("g1Tilde", t, x, k, ans);
	
	return ans;
}

lli comp_g1Tilde(int t, int x, int k) {
	//base case k=0 return 0
	if(k == 0) return 0;
	
	// nonzero k, t=0, also 0
	if(t == 0) return 0;
	
	// nonzero t and k, x zero. Shoulf be undefined?? 
	if(x == 0) {
		// EXCEPTION
		exit(1);
	}
	
	lli ans = 0;
	
	// guess the label set of L
	for(int l = 1; l <= k; ++l) {
		ans += choose(k, l) * f(t, x, l, k - l);
	}
	
	return ans;
}

// at least two comp of G-X, evap at exactly t. Sees all of X.
// Disallowing x=0 
lli g2Tilde(int t, int x, int k) {
		
		if(verbatim) debugIn("g2Tilde", t, x, k);
		
		rec_depth++;
		
		if(table_g2Tilde[t][x][k] == notSaved) table_g2Tilde[t][x][k] = comp_g2Tilde(t, x, k);
		lli ans = table_g2Tilde[t][x][k];
		
		rec_depth--;
		
		if(verbatim) debugOut("g2Tilde", t, x, k, ans);
		
		return ans;
}

lli comp_g2Tilde(int t, int x, int k) {
	//base case k=0 return 0
	if(k == 0) return 0;
	
	// nonzero k, t=0, also 0)
	if(t == 0) return 0;
	
	// nonzero t and k, x zero. Shoulf be undefined?? 
	if(x == 0) {
		// EXCEPTION
		exit(1);
	}
	
	lli ans = 0;
	
	// guess the label set of comp of G-X with lowest label
	for(int kk = 1; kk < k; ++kk) {
		ans += choose(k - 1, kk - 1) * g1Tilde(t, x, kk) * ( g1Tilde(t, x, k-kk) + g2Tilde(t, x, k-kk) );
	}
	
	return ans;
}

// Connected chordal graphs evap at exactly t, |X|=x, L=l, |rest|=k
// Require G-X connected and L union X a clique. 
// allow x = 0
// empty graph (or empty nonX -- i.e l+k=0) evaps at time t=0, otherwise l must be non zero. If t==1 then k must be 0, for t>= 2 k must be non zero.
lli f(int t, int x, int l, int k) {
		
		if(verbatim) debugIn("f", t, x, l, k);
		
		rec_depth++;
		
		if(table_f[t][x][l][k] == notSaved) table_f[t][x][l][k] = comp_f(t, x, l, k);
		lli ans = table_f[t][x][l][k];
		
		rec_depth--;
		
		if(verbatim) debugOut("f", t, x, l, k, ans);
		
		return ans;
}

lli comp_f(int t, int x, int l, int k) {
	// If size of clique is too large, no dice.
	if(x + l > w) return 0;

	// base case, t=0
	if(t == 0) return 0;
	// if(t == 0) return l+k==0 ? 1 : 0;
	
	// additional base case, t=1. Then L must be non zero and k must be zero. x can be anything, including 0 
	if(t == 1) return k == 0 ? 1 : 0;
	
	// t >= 2 then k must be non zero.
	if(k == 0) return 0;
	
	// t >= 2, k >= 1, l >= 1
	lli ans = 0;
	
	// guess the set of labels in components that evaporate at time exactly t-1. 
	// this set must be non zero.
	for(int kk = 1; kk<=k; ++kk) {
		ans += choose(k, kk) * fTilde(t, x, l, kk) * g(t-2, x+l, x, k-kk);
	}
	
	return ans;
}

// Same as f but now every component of G-(X u L) evaporates at time exactly t-1
// Require at least one such component!
// Require G-X connected and L union X a clique. 
// allow x = 0
lli fTilde(int t, int x, int l, int k) {

	if(verbatim) debugIn("fTilde", t, x, l, k);

	rec_depth++;
	
	if(table_fTilde[t][x][l][k] == notSaved) table_fTilde[t][x][l][k] = comp_fTilde(t, x, l, k);
	lli ans = table_fTilde[t][x][l][k];
	
	rec_depth--;

	if(verbatim) debugOut("fTilde", t, x, l, k, ans);

	return ans;
}

lli comp_fTilde(int t, int x, int l, int k) {
	// base case, t=0.
	if(t == 0) return 0;
	
	// second base case t=1, then comps of G-L union X evap at time 0, this is not possible.
	if(t == 1) return 0;
	
	// since t>= 2, k must be non zero.
	if(k == 0) return 0;	
	
	// t >= 2, k >= 1, l >= 1
	lli ans = 0;
	
	ans += fHat(t, x, l, k);
	
	// precisely one comp of G-(X u L) sees all of X u L
	// since l is non zero, at least one vertex is not in the comp, so kk < k
	for(int kk = 1; kk < k; ++kk) {
		ans += choose(k, kk) * g1Tilde(t-1, x+l, kk) * fHat(t,x,l,k-kk);
	}

	// at least 2 comps of G-(X u L) see all of X u L.
	// at least 1 vertex in those componnets. Possibly all vertices. 
	for(int kk = 1; kk <= k; ++kk) {
		ans += choose(k, kk) * g2Tilde(t - 1, x + l, kk) * gHat(t - 1, x + l, x, k - kk);
	}

	return ans;
}

// Same as fTilde but now no component sees all of X union L
// Require at least one such component of G - (X union L)
// Require G-X connected and L union X a clique. 
// allow x = 0
lli fHat(int t, int x, int l, int k) {
	
	if(verbatim) debugIn("fHat", t, x, l, k);
	
	rec_depth++;
	
	if(table_fHat[t][x][l][k] == notSaved) table_fHat[t][x][l][k] = comp_fHat(t, x, l, k);
	lli ans = table_fHat[t][x][l][k];

	rec_depth--;
	
	if(verbatim) debugOut("fHat", t, x, l, k, ans);
	
	return ans;
}

lli comp_fHat(int t, int x, int l, int k) {
	return fHat(t, x, x, l, k);
}

// Overload of fHat but now with z -- every comp of G-X must see someone outside of 1...z
lli fHat(int t, int x, int z, int l, int k) {
	
	if(verbatim) debugIn("fHat", t, x, z, l, k);
	
	rec_depth++;
	
	if(table_fHat5[t][x][z][l][k] == notSaved) table_fHat5[t][x][z][l][k] = comp_fHat(t, x, z, l, k);
	lli ans =  table_fHat5[t][x][z][l][k];
	
	rec_depth--;
	
	if(verbatim) debugOut("fHat", t, x, z, l, k, ans);
	
	return ans;
}

lli comp_fHat(int t, int x, int z, int l, int k) {
	
	// base case, t=0.
	if(t == 0) return 0;
	
	// second base case t=1, then comps of G-L union X evap at time 0, this is not possible. 
	if(t == 1) return 0;
	
	// since t>= 2, k must be non zero.
	if(k == 0) return 0;
	
	// t >= 2, k >= 1, l >= 1
	lli ans = 0;
	
	// guess component comtaining the lowest labeled vertex
	for(int kk = 1; kk <= k; ++kk) {
		for(int xx = 0; xx <= x; ++xx) {
			for(int ll = 0; ll <= l; ++ll) {
				if(xx+ll == 0 || xx+ll == x+l) continue;
				
				lli prod = 1;
				prod *= choose(k-1, kk-1) * choose(l, ll);
				// prod *= (ll > 0 || z == 0) ? choose(x, xx) : choose(x, xx) - choose(z, xx);
				prod *= (ll > 0) ? choose(x, xx) : choose(x, xx) - choose(z, xx);
				prod *= g1Tilde(t-1, xx + ll, kk);
				prod *= ll < l ? fHat(t, x+ll, z, l-ll, k-kk) : gHat(t-1, x+ll, z, k-kk );
				
				ans += prod;
				
			}
		}
	}

	return ans;

}

int main() {
	int n;
	cout << "First input is the number of vertices, second is the bound on the clique size." << endl;
	cin >> n >> w;
	
	chooseInit(n);

	table_chordal.resize(n+1, notSaved);
	table_chordal_conn.resize(n+1, notSaved);
	table_g.resize(n+1, vector<vector<vector<lli> > > (n+1, vector<vector<lli> > (n+1, vector<lli> (n+1, notSaved))));
	table_gTilde.resize(n+1, vector<vector<vector<lli> > > (n+1, vector<vector<lli> > (n+1, vector<lli> (n+1, notSaved))));
	table_gHat.resize(n+1, vector<vector<vector<lli> > > (n+1, vector<vector<lli> > (n+1, vector<lli> (n+1, notSaved))));
	table_g1Tilde.resize(n+1, vector<vector<lli> > (n+1, vector<lli> (n+1, notSaved)));
	table_g2Tilde.resize(n+1, vector<vector<lli> > (n+1, vector<lli> (n+1, notSaved)));
	table_f.resize(n+1, vector<vector<vector<lli> > > (n+1, vector<vector<lli> > (n+1, vector<lli> (n+1, notSaved))));
	table_fTilde.resize(n+1, vector<vector<vector<lli> > > (n+1, vector<vector<lli> > (n+1, vector<lli> (n+1, notSaved))));
	table_fHat.resize(n+1, vector<vector<vector<lli> > > (n+1, vector<vector<lli> > (n+1, vector<lli> (n+1, notSaved))));
	table_fHat5.resize(n+1, vector<vector<vector<vector <lli> > > > (n+1, vector<vector<vector<lli> > > (n+1, vector<vector<lli> > (n+1, vector<lli> (n+1, notSaved)))));
	
	cout << "number of connected labeled chordal graphs: " << chordal_conn(n) << endl;
	cout << "number of labeled chordal graphs: " << chordal(n) << endl;
}
back to top