Commit 25e97342 authored by Szabolcs Horvát's avatar Szabolcs Horvát
Browse files

Initial commit

parent 1883e895
(* Specify build settings such as compiler and linker flags, libraries to be linked, etc. here *)
Switch[$OperatingSystem,
"MacOSX", (* Compilation settings for OS X *)
$buildSettings = {
"CompileOptions" -> {"-std=c++14"}
(*
, "IncludeDirectories" -> {}
, "LibraryDirectories" -> {}
*)
},
"Unix", (* Compilation settings for Linux *)
$buildSettings = {
"CompileOptions" -> {"-std=c++14"}
(*
, "IncludeDirectories" -> {}
, "LibraryDirectories" -> {}
*)
},
"Windows", (* Compilation settings for Windows *)
$buildSettings = {
"CompileOptions" -> {"/EHsc", "/wd4244", "/DNOMINMAX"}
(*
, "IncludeDirectories" -> {}
, "LibraryDirectories" -> {}
*)
}
]
(* Mathematica Package *)
(* Based on LTemplate, https://github.com/szhorvat/LTemplate *)
(* :Title: ConnectedGraphSampler *)
(* :Context: ConnectedGraphSampler` *)
(* :Author: Szabolcs Horvát *)
(* :Date: 2020-08-24 *)
(* :Package Version: 0.1 *)
(* :Mathematica Version: 10.0 *)
(* :Copyright: (c) 2020 Szabolcs Horvát *)
(* :Keywords: *)
(* :Discussion: *)
BeginPackage["ConnectedGraphSampler`"]
(* Privately load LTemplate. Note the leading ` character which ensures that the embedded LTemplate gets loaded. *)
Get["LTemplate`LTemplatePrivate`"];
(* ConfigureLTemplate[] *must* be called at this point. *)
ConfigureLTemplate[
"MessageSymbol" -> ConnectedGraphSampler,
(* If lazy loading is enabled, functions are loaded only on first use.
This improves package loading performance, but it is not convenient
during development and debugging. *)
"LazyLoading" -> False
];
(* Public package symbols go here: *)
ConnectedGraphSampler::usage = "ConnectedGraphSampler is a symbol to which package messages are associated.";
CGSSample::usage =
"CGSSample[degrees] generates a biased sample from the set of graphs with the given degrees. The result is returned in the form {graph, Log[samplingWeight]}.\n" <>
"CGSSample[degrees, Connected -> True] samples only connected graphs.\n" <>
"CGSSample[degrees, n] generates n biased samples.";
CGSSampleProp::usage =
"CGSSampleProp[degrees, prop, n] generates n random graphs of the given degrees, computed prop[graph] for each, then returns the obtained weighted samples as a WeightedData. Use Connected -> True to sample only connected graphs.";
CGSGraphicalQ::usage = "CGSGraphicalQ[degrees] tests if degrees are realized by a simple graph.";
CGSPotentiallyConnectedQ::usage = "CGSPotentiallyConnectedQ[degrees] tests if degrees are potentially connected.";
Connected
`Developer`Recompile::usage = "ConnectedGraphSampler`Developer`Recompile[] recompiles the ConnectedGraphSampler library and reloads the functions.";
PrependTo[$ContextPath, $Context <> "Developer`"];
Begin["`Private`"]
(* Helper function: abort package loading and leave a clean $ContextPath behind *)
packageAbort[] := (End[]; EndPackage[]; Abort[])
$packageVersion = "0.1";
$packageDirectory = DirectoryName[$InputFileName];
$systemID = $SystemID;
(* On OS X libraries use libc++ ABI since M10.4 and libstdc++ ABI up to M10.3.
We need separate binaries to support M10.3 and earlier.
http://community.wolfram.com/groups/-/m/t/816033 *)
If[$OperatingSystem === "MacOSX", $systemID = $systemID <> If[$VersionNumber <= 10.3, "-libstdc++", "-libc++"]];
$libraryDirectory = FileNameJoin[{$packageDirectory, "LibraryResources", $systemID}];
$sourceDirectory = FileNameJoin[{$packageDirectory, "LibraryResources", "Source"}];
$buildSettingsFile = FileNameJoin[{$packageDirectory, "BuildSettings.m"}];
(* Add $libraryDirectory to $LibraryPath in case package is not installed in $UserBaseDirectory/Applications. *)
If[Not@MemberQ[$LibraryPath, $libraryDirectory],
PrependTo[$LibraryPath, $libraryDirectory]
]
(***** The library template *****)
template =
LClass["ConnectedGraphSampler",
{
LFun["setDS", {{Integer, 1, "Constant"}}, "Void"],
LFun["getDS", {}, {Integer, 1}],
LFun["seed", {Integer}, "Void"],
LFun["generateSample", {}, {Integer, 2}],
LFun["generateConnSample", {}, {Integer, 2}],
LFun["getEdges", {}, {Integer, 2}],
LFun["getLogProb", {}, Real],
LFun["graphicalQ", {}, True|False]
}
];
(***** Compilation, loading and initialization *****)
$buildSettings = None;
If[FileExistsQ[$buildSettingsFile], Get[$buildSettingsFile] ]
Recompile::build = "No build settings found. Please check BuildSettings.m."
Recompile[] :=
Module[{},
(* abort compilation if no build settings are present *)
If[$buildSettings === None,
Message[Recompile::build];
Return[$Failed]
];
(* create directory for binary if it doesn't exist yet *)
If[Not@DirectoryQ[$libraryDirectory],
CreateDirectory[$libraryDirectory]
];
(* compile code *)
SetDirectory[$sourceDirectory];
CompileTemplate[template, { (* TODO add any extra .cpp source files to be included in the compilation *) },
"ShellCommandFunction" -> Print, "ShellOutputFunction" -> Print,
"TargetDirectory" -> $libraryDirectory,
Sequence @@ $buildSettings
];
ResetDirectory[];
(* load library *)
loadConnectedGraphSampler[] (* defined below *)
]
loadConnectedGraphSampler[] :=
Module[{deps},
(* mechanism for loading shared library dependencies, if necessary *)
deps = FileNameJoin[{$libraryDirectory, "dependencies.m"}];
Check[
If[FileExistsQ[deps], Get[deps]],
Return[$Failed]
];
(* load the library *)
If[LoadTemplate[template] === $Failed,
Return[$Failed]
];
(* TODO add any other post-loading initialization if necessary *)
]
(* Load library, compile if necessary. *)
If[LoadTemplate[template] === $Failed,
Print[Style["Loading failed, trying to recompile ...", Red]];
If[Recompile[] === $Failed
,
Print[Style["Cannot load or compile library. \[SadSmiley] Aborting.", Red]];
packageAbort[] (* cleanly abort package loading *)
,
Print[Style["Successfully compiled and loaded the library. \[HappySmiley]", Red]];
]
]
(***** Definitions of package functions *****)
cgsTag::usage = "";
throw::usage = "throw[val]";
throw[val_] := Throw[val, cgsTag]
catch::usage = "catch[expr]";
SetAttributes[catch, HoldFirst]
catch[expr_] := Catch[expr, cgsTag]
check::usage = "check[val]";
check[val_LibraryFunctionError] := throw[$Failed] (* this was originally throw[val] *)
check[$Failed] := throw[$Failed]
check[HoldPattern[LibraryFunction[___][___]]] := throw[$Failed]
check[val_] := val
toGraph[n_, opt : OptionsPattern[]][edges_] := Graph[Range[n], edges + 1, Sequence@@FilterRules[{opt}, Options[Graph]]]
CGSGraphicalQ[degrees : {___Integer}] :=
catch@Block[{sampler = Make["ConnectedGraphSampler"]},
If[Max[degrees] >= Length[degrees],
False,
sampler@"setDS"[degrees];
sampler@"graphicalQ"[]
]
]
CGSPotentiallyConnectedQ[degrees : {___Integer}] :=
Total[degrees]/2 >= Length[degrees] - 1
Options[CGSSample] =
Join[
Options[Graph],
{ Connected -> False }
];
SyntaxInformation[CGSSample] = {"ArgumentsPattern" -> {_, _., OptionsPattern[]}};
CGSSample[degrees_, n_Integer ? NonNegative, opt : OptionsPattern[]] :=
catch@Block[{sampler = Make["ConnectedGraphSampler"]},
check@sampler@"setDS"[degrees];
If[OptionValue[Connected],
Table[
{toGraph[Length[degrees], opt]@check@sampler@"generateConnSample"[], sampler@"getLogProb"[]},
{n}
]
,
Table[
{toGraph[Length[degrees], opt]@check@sampler@"generateSample"[], sampler@"getLogProb"[]},
{n}
]
]
]
CGSSample[degrees_, opt : OptionsPattern[]] := catch@First@check@CGSSample[degrees, 1, opt]
Options[CGSSampleProp] = {
Connected -> False
};
SyntaxInformation[CGSSampleProp] = {"ArgumentsPattern" -> {_, _, _, OptionsPattern[]}};
CGSSampleProp[degrees_, prop_, n_Integer ? NonNegative, opt : OptionsPattern[]] :=
catch@Module[{sampler = Make["ConnectedGraphSampler"], values, weights},
check@sampler@"setDS"[degrees];
{values, weights} = Transpose@If[OptionValue[Connected],
Table[
{prop@toGraph[Length[degrees]]@check@sampler@"generateConnSample"[], sampler@"getLogProb"[]},
{n}
]
,
Table[
{prop@toGraph[Length[degrees]]@check@sampler@"generateSample"[], sampler@"getLogProb"[]},
{n}
]
];
weights = weights - Min[weights];
WeightedData[values, Exp[-weights]]
]
End[] (* `Private` *)
EndPackage[]
If[Not@OrderedQ[{10.0, 2}, {$VersionNumber, $ReleaseNumber}],
Print["ConnectedGraphSampler requires Mathematica 10.0.2 or later. Aborting."];
Abort[]
]
Unprotect["ConnectedGraphSampler`*", "ConnectedGraphSampler`Developer`*"];
Get["ConnectedGraphSampler`ConnectedGraphSampler`"]
(* Protect all package symbols *)
SetAttributes[
Evaluate@Flatten[Names /@ {"ConnectedGraphSampler`*", "ConnectedGraphSampler`Developer`*"}],
{Protected, ReadProtected}
]
#include <LTemplate.h>
#define Assert massert
#include "../../../../src/Sampler.h"
#include "../../../../src/ConnSampler.h"
#include <random>
using namespace CDS;
class ConnectedGraphSampler {
DegreeSequence *ds;
std::mt19937 rng;
edgelist_t edges;
double logprob;
public:
ConnectedGraphSampler() :
rng{std::random_device{}()},
logprob(0),
ds(new DegreeSequence)
{ }
~ConnectedGraphSampler() { delete ds; }
void seed(mint s) { rng.seed(s); }
void setDS(mma::IntTensorRef degseq) {
auto old_ds = ds;
ds = new DegreeSequence(degseq.begin(), degseq.end());
delete old_ds;
edges.clear();
logprob = 0;
}
mma::IntTensorRef getDS() const {
return mma::makeVector<mint>(ds->degrees().size(), ds->degrees().data());
}
mma::IntTensorRef getDistrib() const {
return mma::makeVector<mint>(ds->degree_distribution().size(), ds->degree_distribution().data());
}
bool graphicalQ() const { return ds->is_graphical(); }
/*
void decrement(int u) { ds->decrement(u); }
void increment(int u) { ds->increment(u); }
void connect(int u, int v) { ds->connect(u, v); }
mint watershed() const { return ds->watershed(); }
*/
auto getEdges() const {
auto res = mma::makeMatrix<mint>(edges.size(), 2);
for (int i=0; i < edges.size(); ++i) {
res(i,0) = edges[i].first;
res(i,1) = edges[i].second;
}
return res;
}
double getLogProb() const {
return logprob;
}
mma::IntMatrixRef generateSample() {
std::tie(edges, logprob) = CDS::sample(*ds, rng);
return getEdges();
}
mma::IntMatrixRef generateConnSample() {
std::tie(edges, logprob) = CDS::sample_conn(*ds, rng);
return getEdges();
}
};
This diff is collapsed.
#ifndef CDS_COMMON_H
#define CDS_COMMON_H
#include <vector>
#ifndef Assert
#include <cassert>
#define Assert assert
#endif
namespace CDS {
typedef int deg_t;
typedef std::pair<int, int> edge;
typedef std::vector<edge> edgelist_t;
template<typename T>
T sqr(const T &x) { return x*x; }
} // namespace CDS
#endif // CDS_COMMON_H
#ifndef CDS_CONN_SAMPLER_H
#define CDS_CONN_SAMPLER_H
#include "Common.h"
#include "DegreeSequence.h"
#include "EquivClass.h"
#include <vector>
#include <stdexcept>
#include <tuple>
#include <random>
namespace CDS {
template<typename RNG>
std::tuple<edgelist_t, double> sample_conn(DegreeSequence ds, RNG &rng) {
// Null graph considered non-connected
if (ds.n == 0)
throw std::invalid_argument("The degree sequence is not potentially connected.");
if (! ds.is_graphical())
throw std::invalid_argument("The degree sequence is not graphical.");
typedef vector<char> bitmask_t;
EquivClass conn_tracker(ds); // connectivity tracker
if (! conn_tracker.is_potentially_connected())
throw std::invalid_argument("The degree sequence is not potentially connected.");
edgelist_t edges;
double logprob = 0;
bitmask_t exclusion(ds.n);
int vertex = 0;
vector<int> allowed;
vector<double> weights;
while (true) {
if (ds[vertex] == 0) { // No more stubs left on current vertex
if (vertex == ds.n - 1) // All vertices have been processed
return {edges, logprob};
// Advance to next vertex and clear exclusion
vertex += 1;
std::fill(exclusion.begin(), exclusion.end(), 0);
continue;
}
allowed.clear();
weights.clear();
// Construct allowed set
{
// Temporarily connect all but one stub of 'vertex' to highest-degree
// non-excluded vertices. All of these are allowed connections.
DegreeSequence work = ds;
int d = ds[vertex];
auto supernode = conn_tracker.get_class(vertex);
int d_supernode = supernode->degree();
int comp_count = conn_tracker.component_count();
int edge_count = conn_tracker.edge_count();
int i=ds.n-1;
while (d > 1) {
int v = ds.sorted_verts[i--];
Assert(work[v] > 0);
if (v != vertex && ! exclusion[v]) {
// mma::mout << "TC: " << vertex << "-" << v << " ";
work.connect(vertex, v);
if ( comp_count == 1 ||
edge_count == 1 ||
(d_supernode > 2 && edge_count > comp_count - 1) ||
(conn_tracker.get_class(v) != supernode && (d_supernode > 1 || conn_tracker.get_class(v)->degree() > 1)) )
{
allowed.push_back(v);
weights.push_back(ds[v]);
}
d--;
}
}
// Remove the final stub of 'vertex'.
work.decrement(vertex);
// Find watershed degree.
int wd = work.watershed();
// Of the rest of the vertices, determine if a connection is allowed
// based on the watershed degree.
for (; i >= 0; --i) {
int v = ds.sorted_verts[i];
if (ds[v] >= wd) {
if (v != vertex && ! exclusion[v]) {
if ( comp_count == 1 ||
edge_count == 1 ||
(d_supernode > 2 && edge_count > comp_count - 1) ||
(conn_tracker.get_class(v) != supernode && (d_supernode > 1 || conn_tracker.get_class(v)->degree() > 1)) )
{
allowed.push_back(v);
weights.push_back(ds[v]);
}
}
} else {
break;
}
}
}
Assert(! allowed.empty());
double tot = std::accumulate(weights.begin(), weights.end(), 0.0);
logprob -= std::log(tot);
std::discrete_distribution<> choose(weights.begin(), weights.end());
int u = allowed[choose(rng)];
exclusion[u] = 1;
ds.connect(u, vertex);
conn_tracker.connect(u, vertex);
edges.push_back({vertex, u});
}
}
} // namespace CDS
#endif // CDS_CONN_SAMPLER_H
#ifndef CDS_DEGREE_SEQUENCE_H
#define CDS_DEGREE_SEQUENCE_H
#include "Common.h"
#include <vector>
#include <numeric>
#include <algorithm>
#include <stdexcept>
namespace CDS {
using std::vector;
class DegreeSequence {
vector<deg_t> degseq; // the degree sequence
const int n; // number of degrees
vector<int> deg_counts; // deg_counts[d] is the number of vertices with degree d
vector<int> accum_counts; // accum_counts[d] is the number of vertices with degree <= d
vector<int> sorted_verts; // vertex indices sorted by vertex degree
vector<int> sorted_index; // sorted_index[u] is the index of vertex u in sorted_verts, i.e. sorted_verts[sorted_index[u]] == u
deg_t dmax, dmin; // largest and smallest non-zero (!) degree, assuming n_nonzero != 0
int n_nonzero; // number of non-zero degrees
int dsum; // the sum of degrees
// d(i) returns d_i in the non-increasingly sorted degree sequence
// note that i is assumed to use 1-based indexing!
deg_t d(int i) const { return degseq[sorted_verts[n - i]]; }
public:
DegreeSequence() : n(0), dmax(0), dmin(0), n_nonzero(0), dsum(0) { }
// Initialize degree sequence, O(n)
template<typename It>
DegreeSequence(It first, It last) :
degseq(first, last),
n(degseq.size()),
deg_counts(n),
sorted_verts(n), sorted_index(n),
accum_counts(n)
{
// initialize sorted_verts
std::iota(sorted_verts.begin(), sorted_verts.end(), 0);
std::sort(sorted_verts.begin(), sorted_verts.end(), [this] (int u, int v) { return degseq[u] < degseq[v]; } );
// initialize sorted_index
for (int i=0; i < n; ++i)
sorted_index[sorted_verts[i]] = i;
dmax = 0;
dmin = 0;
n_nonzero = 0;
dsum = 0;
for (int i=0; i < n; ++i) {
deg_t d = degseq[i];
if (d < 0)
throw std::domain_error("Degrees must be non-negative.");