Commit 28bf0417 authored by Szabolcs Horvát's avatar Szabolcs Horvát
Browse files

Add multigraph sampler and cleanup

parent 7069d293
......@@ -4,7 +4,7 @@
Switch[$OperatingSystem,
"MacOSX", (* Compilation settings for OS X *)
$buildSettings = {
"CompileOptions" -> {"-std=c++14", "-O3", "-march=native"}
"CompileOptions" -> {"-std=c++14", "-O3"}
(*
, "IncludeDirectories" -> {}
......@@ -14,7 +14,7 @@ Switch[$OperatingSystem,
"Unix", (* Compilation settings for Linux *)
$buildSettings = {
"CompileOptions" -> {"-std=c++14", "-O3", "-march=native"}
"CompileOptions" -> {"-std=c++14", "-O3"}
(*
, "IncludeDirectories" -> {}
......
......@@ -32,16 +32,23 @@ ConnectedGraphSampler::usage = "ConnectedGraphSampler is a symbol to which packa
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.";
CGSSampleWeights::usage = "CGSSampleWeights[degrees, n]";
"CGSSample[degrees, n] generates n biased samples.\n" <>
"CGSSample[degrees, \"Connected\" -> True] samples only connected graphs.\n" <>
"CGSSample[degrees, \"MultiEdges\" -> True] allows multi-edges.\n" <>
"CGSSample[degrees, Exponent -> \[Alpha]] sets the degree affinity exponent.";
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.";
CGSSamplePropRaw::usage = "CGSSamplePropRaw[degrees, prop, n] generates n random graphs of the given degrees, computed prop[graph] for each, then returns the obtained weighted samples as a list of {value, Log[samplingWeight]} pairs. Use Connected -> True to sample only connected graphs.";
CGSToWeightedData::usage = "CGSToWeightedData[samples] converts a list of {value, Log[samplingWeight]} pairs to a WeightedData expression.";
"CGSSampleProp[degrees, prop, n] generates n random graphs with the given degrees, computes prop[graph] for each, then returns the obtained weighted samples as a WeightedData. Accepts the same options as CGSSample.";
CGSSampleWeights::usage = "CGSSampleWeights[degrees, n] generates n random graphs with the given degrees and returns the logarithms of their sampling weights. Accepts the same options as CGSSample.";
CGSSamplePropRaw::usage = "CGSSamplePropRaw[degrees, prop, n] generates n random graphs with the given degrees, computes value = prop[graph] for each, and returns the result as {value, Log[samplingWeight]} pairs. Accepts the same options as CGSSample.";
CGSToWeightedData::usage = "CGSToWeightedData[rawData] converts a list of {value, Log[samplingWeight]} pairs to a WeightedData expression.";
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.";
......@@ -76,16 +83,32 @@ If[Not@MemberQ[$LibraryPath, $libraryDirectory],
(***** The library template *****)
template =
LClass["ConnectedGraphSampler",
LTemplate[
"ConnectedGraphSampler",
{
LFun["setDS", {{Integer, 1, "Constant"} (* degree sequence *)}, "Void"],
LFun["getDS", {}, {Integer, 1}],
LFun["seed", {Integer}, "Void"],
LFun["generateSample", {Real (* alpha *)}, {Integer, 2}],
LFun["generateConnSample", {Real (* alpha *)}, {Integer, 2}],
LFun["getEdges", {}, {Integer, 2}],
LFun["getLogProb", {}, Real],
LFun["graphicalQ", {}, True|False]
LClass["ConnectedGraphSampler",
{
LFun["setDS", {{Integer, 1, "Constant"} (* degree sequence *)}, "Void"],
LFun["getDS", {}, {Integer, 1}],
LFun["seed", {Integer}, "Void"],
LFun["generateSample", {Real (* alpha *)}, {Integer, 2}],
LFun["generateConnSample", {Real (* alpha *)}, {Integer, 2}],
LFun["getEdges", {}, {Integer, 2}],
LFun["getLogProb", {}, Real],
LFun["graphicalQ", {}, True | False]
}
],
LClass["ConnectedGraphSamplerMulti",
{
LFun["setDS", {{Integer, 1, "Constant"} (* degree sequence *)}, "Void"],
LFun["getDS", {}, {Integer, 1}],
LFun["seed", {Integer}, "Void"],
LFun["generateSample", {Real (* alpha *)}, {Integer, 2}],
LFun["generateConnSample", {Real (* alpha *)}, {Integer, 2}],
LFun["getEdges", {}, {Integer, 2}],
LFun["getLogProb", {}, Real]
}
]
}
];
......@@ -151,9 +174,9 @@ If[LoadTemplate[template] === $Failed,
]
(***** Definitions of package functions *****)
(***** Helper functions *****)
cgsTag::usage = "";
cgsTag::usage = "Tag used by throw/catch.";
throw::usage = "throw[val]";
throw[val_] := Throw[val, cgsTag]
......@@ -163,13 +186,16 @@ SetAttributes[catch, HoldFirst]
catch[expr_] := Catch[expr, cgsTag]
check::usage = "check[val]";
check[val_LibraryFunctionError] := throw[$Failed] (* this was originally throw[val] *)
check[val_LibraryFunctionError] := throw[$Failed]
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]]]
(***** Definitions of package functions *****)
CGSGraphicalQ[degrees : {___Integer}] :=
catch@Block[{sampler = Make["ConnectedGraphSampler"]},
If[Max[degrees] >= Length[degrees],
......@@ -179,24 +205,29 @@ CGSGraphicalQ[degrees : {___Integer}] :=
]
]
CGSPotentiallyConnectedQ[degrees : {___Integer}] :=
Total[degrees]/2 >= Length[degrees] - 1
CGSPotentiallyConnectedQ[{0}] := True (* special case of a single isolated vertex *)
CGSPotentiallyConnectedQ[degrees : {___Integer ? Positive}] :=
Total[degrees]/2 >= Length[degrees] - 1 && FreeQ[degrees, 0, {1}]
CGSPotentiallyConnectedQ[_] := False
Options[CGSSample] =
Join[
Options[Graph],
{
Connected -> False,
"MultiEdges" -> False,
"Connected" -> False,
RandomSeeding -> Automatic,
Exponent -> 1
}
];
SyntaxInformation[CGSSample] = {"ArgumentsPattern" -> {_, _., OptionsPattern[]}};
CGSSample[degrees_, n_Integer ? NonNegative, opt : OptionsPattern[]] :=
catch@Block[{sampler = Make["ConnectedGraphSampler"]},
catch@Block[{sampler = If[TrueQ@OptionValue["MultiEdges"], Make["ConnectedGraphSamplerMulti"], Make["ConnectedGraphSampler"]]},
check@sampler@"setDS"[degrees];
sampler@"seed"[ Replace[OptionValue[RandomSeeding], Automatic :> RandomInteger[2^31-1]] ];
If[OptionValue[Connected],
If[TrueQ@OptionValue["Connected"],
Table[
{toGraph[Length[degrees], opt]@check@sampler@"generateConnSample"[OptionValue[Exponent]], sampler@"getLogProb"[]},
{n}
......@@ -212,15 +243,16 @@ CGSSample[degrees_, opt : OptionsPattern[]] := catch@First@check@CGSSample[degre
Options[CGSSampleWeights] = {
Connected -> False,
"MultiEdges" -> False,
"Connected" -> False,
RandomSeeding -> Automatic,
Exponent -> 1
};
CGSSampleWeights[degrees_, n_Integer ? NonNegative, opt : OptionsPattern[]] :=
catch@Block[{sampler = Make["ConnectedGraphSampler"]},
catch@Block[{sampler = If[TrueQ@OptionValue["MultiEdges"], Make["ConnectedGraphSamplerMulti"], Make["ConnectedGraphSampler"]]},
check@sampler@"setDS"[degrees];
sampler@"seed"[ Replace[OptionValue[RandomSeeding], Automatic :> RandomInteger[2^31-1]] ];
If[OptionValue[Connected],
If[TrueQ@OptionValue["Connected"],
Table[
check@sampler@"generateConnSample"[OptionValue[Exponent]];
sampler@"getLogProb"[],
......@@ -237,7 +269,8 @@ CGSSampleWeights[degrees_, n_Integer ? NonNegative, opt : OptionsPattern[]] :=
Options[CGSSampleProp] = {
Connected -> False,
"MultiEdges" -> False,
"Connected" -> False,
RandomSeeding -> Automatic,
Exponent -> 1
};
......@@ -247,16 +280,17 @@ CGSSampleProp[degrees_, prop_, n_Integer ? NonNegative, opt : OptionsPattern[]]
Options[CGSSamplePropRaw] = {
Connected -> False,
"MultiEdges" -> False,
"Connected" -> False,
RandomSeeding -> Automatic,
Exponent -> 1
};
SyntaxInformation[CGSSamplePropRaw] = {"ArgumentsPattern" -> {_, _, _, OptionsPattern[]}};
CGSSamplePropRaw[degrees_, prop_, n_Integer ? NonNegative, opt : OptionsPattern[]] :=
catch@Block[{sampler = Make["ConnectedGraphSampler"]},
catch@Block[{sampler = If[TrueQ@OptionValue["MultiEdges"], Make["ConnectedGraphSamplerMulti"], Make["ConnectedGraphSampler"]]},
check@sampler@"setDS"[degrees];
sampler@"seed"[ Replace[OptionValue[RandomSeeding], Automatic :> RandomInteger[2^31-1]] ];
If[OptionValue[Connected],
If[TrueQ@OptionValue["Connected"],
Table[
{prop@toGraph[Length[degrees]]@check@sampler@"generateConnSample"[OptionValue[Exponent]], sampler@"getLogProb"[]},
{n}
......@@ -269,6 +303,7 @@ CGSSamplePropRaw[degrees_, prop_, n_Integer ? NonNegative, opt : OptionsPattern[
]
]
CGSToWeightedData[data_] :=
Module[{weights, values},
{values, weights} = Transpose[data];
......@@ -276,6 +311,7 @@ CGSToWeightedData[data_] :=
WeightedData[values, Exp[-weights]]
]
End[] (* `Private` *)
EndPackage[]
#ifndef CONNECTED_GRAPH_SAMPLER
#define CONNECTED_GRAPH_SAMPLER
#include <LTemplate.h>
#define Assert massert
......@@ -12,8 +15,8 @@ using namespace CDS;
class ConnectedGraphSampler {
DegreeSequence *ds;
std::mt19937 rng;
DegreeSequence *ds;
edgelist_t edges;
double logprob;
......@@ -22,8 +25,8 @@ public:
ConnectedGraphSampler() :
rng{std::random_device{}()},
logprob(0),
ds(new DegreeSequence)
ds(new DegreeSequence),
logprob(0)
{ }
~ConnectedGraphSampler() { delete ds; }
......@@ -80,3 +83,5 @@ public:
return getEdges();
}
};
#endif // CONNECTED_GRAPH_SAMPLER
#ifndef CONNECTED_GRAPH_SAMPLER_MULTI
#define CONNECTED_GRAPH_SAMPLER_MULTI
#include <LTemplate.h>
#define Assert massert
#include "../../../../src/MultiSampler.h"
#include "../../../../src/ConnMultiSampler.h"
#include <random>
using namespace CDS;
class ConnectedGraphSamplerMulti {
std::mt19937 rng;
std::vector<deg_t> ds;
edgelist_t edges;
double logprob;
public:
ConnectedGraphSamplerMulti() :
rng{std::random_device{}()},
logprob(0)
{ }
void seed(mint s) { rng.seed(s); }
void setDS(mma::IntTensorRef degseq) {
ds.assign(degseq.begin(), degseq.end());
edges.clear();
logprob = 0;
}
mma::IntTensorRef getDS() const {
return mma::makeVector<mint>(ds.size(), ds.data());
}
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(double alpha) {
std::tie(edges, logprob) = CDS::sample_multi(ds, alpha, rng);
return getEdges();
}
mma::IntMatrixRef generateConnSample(double alpha) {
std::tie(edges, logprob) = CDS::sample_conn_multi(ds, alpha, rng);
return getEdges();
}
};
#endif // CONNECTED_GRAPH_SAMPLER_MULTI
This diff is collapsed.
......@@ -2,6 +2,7 @@
#define CDS_COMMON_H
#include <vector>
#include <cmath>
#ifndef Assert
#include <cassert>
......@@ -17,9 +18,282 @@ typedef std::vector<edge> edgelist_t;
typedef std::vector<char> bitmask_t;
template<typename T>
T sqr(const T &x) { return x*x; }
// Compute the logarithm of n-factorial.
// Based on: https://www.johndcook.com/blog/2010/08/16/how-to-compute-log-factorial/
inline double logfact(int n) {
const double lf[] = {
0.0,
0.0,
0.6931471805599453,
1.791759469228055,
3.1780538303479458,
4.787491742782046,
6.579251212010101,
8.525161361065415,
1.060460290274525e1,
1.2801827480081469e1,
1.5104412573075516e1,
1.7502307845873887e1,
1.9987214495661885e1,
2.2552163853123425e1,
2.519122118273868e1,
2.789927138384089e1,
3.0671860106080672e1,
3.350507345013689e1,
3.639544520803305e1,
3.9339884187199495e1,
4.2335616460753485e1,
4.538013889847691e1,
4.847118135183523e1,
5.160667556776438e1,
5.478472939811232e1,
5.800360522298052e1,
6.1261701761002e1,
6.455753862700634e1,
6.788974313718154e1,
7.125703896716801e1,
7.465823634883016e1,
7.80922235533153e1,
8.155795945611504e1,
8.505446701758152e1,
8.858082754219768e1,
9.21361756036871e1,
9.57196945421432e1,
9.933061245478743e1,
1.0296819861451381e2,
1.0663176026064346e2,
1.1032063971475739e2,
1.140342117814617e2,
1.1777188139974507e2,
1.2153308151543864e2,
1.253172711493569e2,
1.2912393363912722e2,
1.3295257503561632e2,
1.3680272263732635e2,
1.4067392364823425e2,
1.445657439463449e2,
1.4847776695177302e2,
1.5240959258449735e2,
1.563608363030788e2,
1.603311282166309e2,
1.6432011226319517e2,
1.6832744544842765e2,
1.723527971391628e2,
1.7639584840699735e2,
1.8045629141754378e2,
1.8453382886144948e2,
1.886281734236716e2,
1.927390472878449e2,
1.9686618167289e2,
2.0100931639928152e2,
2.051681994826412e2,
2.0934258675253685e2,
2.1353224149456327e2,
2.1773693411395422e2,
2.2195644181913033e2,
2.261905483237276e2,
2.3043904356577696e2,
2.3470172344281826e2,
2.3897838956183432e2,
2.432688490029827e2,
2.4757291409618688e2,
2.518904022097232e2,
2.5622113555000954e2,
2.605649409718632e2,
2.649216497985528e2,
2.692910976510198e2,
2.736731242856937e2,
2.780675734403661e2,
2.824742926876304e2,
2.86893133295427e2,
2.913239500942703e2,
2.9576660135076065e2,
3.0022094864701415e2,
3.046868567656687e2,
3.091641935801469e2,
3.1365282994987905e2,
3.181526396202093e2,
3.2266349912672615e2,
3.271852877037752e2,
3.317178871969285e2,
3.3626118197919845e2,
3.40815058870799e2,
3.4537940706226686e2,
3.4995411804077025e2,
3.545390855194408e2,
3.591342053695754e2,
3.6373937555556347e2,
3.6835449607240474e2,
3.72979468885689e2,
3.7761419787391867e2,
3.8225858877306e2,
3.8691254912321756e2,
3.915759882173296e2,
3.9624881705179155e2,
4.0093094827891576e2,
4.056222961611449e2,
4.1032277652693733e2,
4.1503230672824964e2,
4.197508055995447e2,
4.244781934182571e2,
4.2921439186665157e2,
4.339593239950148e2,
4.3871291418612117e2,
4.4347508812091894e2,
4.482457727453846e2,
4.530248962384961e2,
4.5781238798127816e2,
4.626081785268749e2,
4.674121995716082e2,
4.722243839269806e2,
4.7704466549258564e2,
4.8187297922988796e2,
4.867092611368394e2,
4.91553448223298e2,
4.9640547848721764e2,
5.012652908915793e2,
5.061328253420349e2,
5.11008022665236e2,
5.158908245878224e2,
5.207811737160441e2,
5.25679013515995e2,
5.305842882944335e2,
5.354969431801695e2,
5.404169241059976e2,
5.453441777911548e2,
5.502786517242855e2,
5.552202941468948e2,
5.60169054037273e2,
5.651248810948744e2,
5.700877257251342e2,
5.750575390247102e2,
5.800342727671308e2,
5.850178793888391e2,
5.900083119756179e2,
5.95005524249382e2,
6.000094705553274e2,
6.050201058494237e2,
6.100373856862386e2,
6.150612662070849e2,
6.200917041284773e2,
6.25128656730891e2,
6.301720818478102e2,
6.352219378550598e2,
6.40278183660408e2,
6.45340778693435e2,
6.504096828956552e2,
6.554848567108891e2,
6.605662610758735e2,
6.65653857411106e2,
6.707476076119127e2,
6.758474740397369e2,
6.809534195136374e2,
6.86065407301994e2,
6.911834011144108e2,
6.96307365093814e2,
7.01437263808737e2,
7.065730622457874e2,
7.1171472580229e2,
7.168622202791034e2,
7.220155118736012e2,
7.271745671728157e2,
7.323393531467393e2,
7.375098371417774e2,
7.426859868743512e2,
7.478677704246434e2,
7.530551562304842e2,
7.582481130813743e2,
7.634466101126401e2,
7.686506167997169e2,
7.738601029525583e2,
7.790750387101673e2,
7.842953945352457e2,
7.895211412089589e2,
7.947522498258135e2,
7.999886917886435e2,
8.05230438803703e2,
8.104774628758636e2,
8.157297363039102e2,
8.209872316759379e2,
8.262499218648428e2,
8.315177800239062e2,
8.367907795824699e2,
8.420688942417004e2,
8.473520979704384e2,
8.52640365001133e2,
8.579336698258575e2,
8.632319871924054e2,
8.685352921004645e2,
8.738435597978657e2,
8.791567657769075e2,
8.844748857707517e2,
8.897978957498901e2,
8.951257719186798e2,
9.004584907119452e2,
9.057960287916464e2,
9.111383630436112e2,
9.164854705743287e2,
9.218373287078048e2,
9.271939149824768e2,
9.325552071481862e2,
9.379211831632081e2,
9.432918211913358e2,
9.486670995990199e2,
9.540469969525603e2,
9.594314920153495e2,
9.648205637451659e2,
9.702141912915183e2,
9.756123539930361e2,
9.810150313749083e2,
9.864222031463685e2,
9.918338491982234e2,
9.97249949600428e2,
1.0026704845997002e3,
1.0080954346171816e3,
1.0135247802461361e3,
1.0189585022496902e3,
1.0243965815586134e3,
1.0298389992691352e3,
1.0352857366408016e3,
1.0407367750943672e3,
1.046192096209725e3,
1.0516516817238692e3,
1.0571155135288948e3,
1.0625835736700299e3,
1.0680558443437014e3,
1.0735323078956328e3,
1.0790129468189748e3,
1.0844977437524656e3,
1.0899866814786221e3,
1.0954797429219627e3,
1.100976911147256e3,
1.1064781693578007e3,
1.111983500893733e3,
1.117492889230361e3,
1.123006317976526e3,
1.1285237708729908e3,
1.134045231790853e3,
1.1395706847299848e3,
1.145100113817496e3,
1.1506335033062237e3,
1.1561708375732421e3,
1.1617121011184006e3,
1.1672572785628802e3