Source code for graph_tool.inference.uncertain_blockmodel
#! /usr/bin/env python# -*- coding: utf-8 -*-## graph_tool -- a general graph manipulation python module## Copyright (C) 2006-2023 Tiago de Paula Peixoto <tiago@skewed.de>## This program is free software; you can redistribute it and/or modify it under# the terms of the GNU Lesser General Public License as published by the Free# Software Foundation; either version 3 of the License, or (at your option) any# later version.## This program 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 Lesser General Public License for more# details.## You should have received a copy of the GNU Lesser General Public License# along with this program. If not, see <http://www.gnu.org/licenses/>.from..import_prop,Graph,GraphView,_get_rng,PropertyMap, \
edge_endpoint_propertyfrom..dl_importimportdl_importdl_import("from . import libgraph_tool_inference as libinference")from.base_statesimport*from.blockmodelimport*from.nested_blockmodelimport*importcollections.abc
[docs]classUncertainBaseState(object):r"""Base state for uncertain network inference."""def__init__(self,g,nested=True,state_args={},bstate=None,self_loops=False,init_empty=False,max_m=1<<16):self.g=gstate_args=dict(state_args)ifbstateisNone:ifinit_empty:self.u=Graph(directed=g.is_directed())self.u.add_vertex(g.num_vertices())self.eweight=self.u.new_ep("int",val=1)elif"g"instate_args:self.u=state_args.pop("g")self.eweight=state_args.pop("eweight",self.u.new_ep("int",val=1))else:self.u=g.copy()self.eweight=self.u.new_ep("int",val=1)else:self.u=bstate.g.copy()ifnested:self.eweight=bstate.levels[0].eweightelse:self.eweight=bstate.eweightself.eweight=self.u.own_property(self.eweight.copy())ifnested:bstate=bstate.copy(g=self.u,state_args=dict(bstate.state_args,eweight=self.eweight))else:bstate=bstate.copy(g=self.u,eweight=self.eweight)self.u.set_fast_edge_removal()self.self_loops=self_loopsN=self.u.num_vertices()ifself.u.is_directed():ifself_loops:M=N*Nelse:M=N*(N-1)else:ifself_loops:M=(N*(N+1))/2else:M=(N*(N-1))/2self.M=MifbstateisNone:ifnested:state_args["state_args"]=state_args.get("state_args",{})state_args["state_args"].update(dict(eweight=self.eweight))self.nbstate=NestedBlockState(self.u,**state_args)self.bstate=self.nbstate.levels[0]else:self.nbstate=Noneself.bstate=BlockState(self.u,eweight=self.eweight,**state_args)else:ifnested:self.nbstate=bstateself.bstate=bstate.levels[0]else:self.nbstate=Noneself.bstate=bstateedges=self.g.get_edges()edges=numpy.concatenate((edges,numpy.ones(edges.shape,dtype=edges.dtype)*(N+1)))self.slist=Vector_size_t(init=edges[:,0])self.tlist=Vector_size_t(init=edges[:,1])self.max_m=max_minit_q_cache()def__getstate__(self):returndict(g=self.g,nested=self.nbstateisnotNone,bstate=(self.nbstateifself.nbstateisnotNoneelseself.bstate),self_loops=self.self_loops,max_m=self.max_m)def__setstate__(self,state):self.__init__(**state)def_get_entropy_args(self,kwargs):ea=self.bstate._get_entropy_args(kwargs,ignore=["latent_edges","density"])uea=libinference.uentropy_args(ea)uea.latent_edges=kwargs.get("latent_edges",True)uea.density=kwargs.get("density",True)returnuea
[docs]defget_block_state(self):"""Return the underlying block state, which can be either :class:`~graph_tool.inference.BlockState` or :class:`~graph_tool.inference.NestedBlockState`. """ifself.nbstateisNone:returnself.bstateelse:returnself.nbstate
[docs]@copy_state_wrapdefentropy(self,latent_edges=True,density=True,**kwargs):"""Return the entropy, i.e. negative log-likelihood."""S=self._state.entropy(latent_edges,density)ifself.nbstateisNone:S+=self.bstate.entropy(**kwargs)else:S+=self.nbstate.entropy(**kwargs)returnS
[docs]@mcmc_sweep_wrapdefmcmc_sweep(self,r=.5,multiflip=True,**kwargs):r"""Perform sweeps of a Metropolis-Hastings acceptance-rejection sampling MCMC to sample network partitions and latent edges. The parameter ``r`` controls the probability with which edge move will be attempted, instead of partition moves. The remaining keyword parameters will be passed to :meth:`~graph_tool.inference.BlockState.mcmc_sweep` or :meth:`~graph_tool.inference.BlockState.multiflip_mcmc_sweep`, if ``multiflip=True``. """ifmultiflip:returnself._algo_sweep(lambdas,**kw:s.multiflip_mcmc_sweep(**kw),r=r,**kwargs)else:returnself._algo_sweep(lambdas,**kw:s.mcmc_sweep(**kw),r=r,**kwargs)
[docs]defmultiflip_mcmc_sweep(self,**kwargs):r"""Alias for :meth:`~UncertainBaseState.mcmc_sweep` with ``multiflip=True``."""returnself.mcmc_sweep(multiflip=True,**kwargs)
[docs]defget_edge_prob(self,u,v,entropy_args={},epsilon=1e-8):r"""Return conditional posterior log-probability of edge :math:`(u,v)`."""ea=self._get_entropy_args(entropy_args)returnself._state.get_edge_prob(u,v,ea,epsilon)
[docs]defget_edges_prob(self,elist,entropy_args={},epsilon=1e-8):r"""Return conditional posterior log-probability of an edge list, with shape :math:`(E,2)`."""ea=self._get_entropy_args(entropy_args)elist=numpy.asarray(elist,dtype="uint64")probs=numpy.zeros(elist.shape[0])self._state.get_edges_prob(elist,probs,ea,epsilon)returnprobs
[docs]defget_graph(self):r"""Return the current inferred graph."""ifself.self_loops:u=GraphView(self.u,efilt=self.eweight.fa>0)else:es=edge_endpoint_property(self.u,self.u.vertex_index,"source")et=edge_endpoint_property(self.u,self.u.vertex_index,"target")u=GraphView(self.u,efilt=numpy.logical_and(self.eweight.fa>0,es.fa!=et.fa))returnu
[docs]defcollect_marginal(self,g=None):r"""Collect marginal inferred network during MCMC runs. Parameters ---------- g : :class:`~graph_tool.Graph` (optional, default: ``None``) Previous marginal graph. Returns ------- g : :class:`~graph_tool.Graph` New marginal graph, with internal edge :class:`~graph_tool.EdgePropertyMap` ``"eprob"``, containing the marginal probabilities for each edge. Notes ----- The posterior marginal probability of an edge :math:`(i,j)` is defined as .. math:: \pi_{ij} = \sum_{\boldsymbol A}A_{ij}P(\boldsymbol A|\boldsymbol D) where :math:`P(\boldsymbol A|\boldsymbol D)` is the posterior probability given the data. """ifgisNone:g=Graph(directed=self.g.is_directed())g.add_vertex(self.g.num_vertices())g.gp.count=g.new_gp("int",0)g.ep.count=g.new_ep("int")if"eprob"noting.ep:g.ep.eprob=g.new_ep("double")u=self.get_graph()libinference.collect_marginal(g._Graph__graph,u._Graph__graph,_prop("e",g,g.ep.count))g.gp.count+=1g.ep.eprob.fa=g.ep.count.fag.ep.eprob.fa/=g.gp.countreturng
[docs]defcollect_marginal_multigraph(self,g=None):r"""Collect marginal latent multigraph during MCMC runs. Parameters ---------- g : :class:`~graph_tool.Graph` (optional, default: ``None``) Previous marginal multigraph. Returns ------- g : :class:`~graph_tool.Graph` New marginal graph, with internal edge :class:`~graph_tool.EdgePropertyMap` ``"w"`` and ``"wcount"``, containing the edge multiplicities and their respective counts. Notes ----- The mean posterior marginal multiplicity distribution of a multi-edge :math:`(i,j)` is defined as .. math:: \pi_{ij}(w) = \sum_{\boldsymbol G}\delta_{w,G_{ij}}P(\boldsymbol G|\boldsymbol D) where :math:`P(\boldsymbol G|\boldsymbol D)` is the posterior probability of a multigraph :math:`\boldsymbol G` given the data. """ifgisNone:g=Graph(directed=self.g.is_directed())g.add_vertex(self.g.num_vertices())g.ep.w=g.new_ep("vector<int>")g.ep.wcount=g.new_ep("vector<int>")libinference.collect_marginal_count(g._Graph__graph,self.u._Graph__graph,_prop("e",self.u,self.eweight),_prop("e",g,g.ep.w),_prop("e",g,g.ep.wcount))returng
[docs]classUncertainBlockState(UncertainBaseState):r"""Inference state of an uncertain graph, using the stochastic block model as a prior. Parameters ---------- g : :class:`~graph_tool.Graph` Measured graph. q : :class:`~graph_tool.EdgePropertyMap` Edge probabilities in range :math:`[0,1]`. q_default : ``float`` (optional, default: ``0.``) Non-edge probability in range :math:`[0,1]`. aE : ``float`` (optional, default: ``NaN``) Expected total number of edges used in prior. If ``NaN``, a flat prior will be used instead. nested : ``boolean`` (optional, default: ``True``) If ``True``, a :class:`~graph_tool.inference.NestedBlockState` will be used, otherwise :class:`~graph_tool.inference.BlockState`. state_args : ``dict`` (optional, default: ``{}``) Arguments to be passed to :class:`~graph_tool.inference.NestedBlockState` or :class:`~graph_tool.inference.BlockState`. bstate : :class:`~graph_tool.inference.NestedBlockState` or :class:`~graph_tool.inference.BlockState` (optional, default: ``None``) If passed, this will be used to initialize the block state directly. self_loops : bool (optional, default: ``False``) If ``True``, it is assumed that the uncertain graph can contain self-loops. References ---------- .. [peixoto-reconstructing-2018] Tiago P. Peixoto, "Reconstructing networks with unknown and heterogeneous errors", Phys. Rev. X 8 041011 (2018). :doi:`10.1103/PhysRevX.8.041011`, :arxiv:`1806.07956` """def__init__(self,g,q,q_default=0.,aE=numpy.nan,nested=True,state_args={},bstate=None,self_loops=False,**kwargs):super().__init__(g,nested=nested,state_args=state_args,bstate=bstate,self_loops=self_loops,**kwargs)self._q=qself._q_default=q_defaultself.p=(q.fa.sum()+(self.M-g.num_edges())*q_default)/self.Mself.q=self.g.new_ep("double",vals=log(q.fa)-log1p(-q.fa))self.q.fa-=log(self.p)-log1p(-self.p)ifq_default>0:self.q_default=log(q_default)-log1p(q_default)self.q_default-=log(self.p)-log1p(-self.p)else:self.q_default=-numpy.infself.S_const=(log1p(-q.fa[q.fa<1]).sum()+log1p(-q_default)*(self.M-self.g.num_edges())-self.M*log1p(-self.p))self.aE=aEifnumpy.isnan(aE):self.E_prior=Falseelse:self.E_prior=Trueself._state=libinference.make_uncertain_state(self.bstate._state,self)def__getstate__(self):state=super().__getstate__()returndict(state,q=self._q,q_default=self._q_default,aE=self.aE)
[docs]defcopy(self,**kwargs):"""Return a copy of the state."""returnUncertainBlockState(**dict(self.__getstate__(),**kwargs))
def__copy__(self):returnself.copy()def__repr__(self):return"<UncertainBlockState object with %s, at 0x%x>"% \
(self.nbstateifself.nbstateisnotNoneelseself.bstate,id(self))def_mcmc_sweep(self,mcmc_state):returnlibinference.mcmc_uncertain_sweep(mcmc_state,self._state,_get_rng())
[docs]classLatentMultigraphBlockState(UncertainBaseState):r"""Inference state of an erased Poisson multigraph, using the stochastic block model as a prior. Parameters ---------- g : :class:`~graph_tool.Graph` Measured graph. aE : ``float`` (optional, default: ``NaN``) Expected total number of edges used in prior. If ``NaN``, a flat prior will be used instead. nested : ``boolean`` (optional, default: ``True``) If ``True``, a :class:`~graph_tool.inference.NestedBlockState` will be used, otherwise :class:`~graph_tool.inference.BlockState`. state_args : ``dict`` (optional, default: ``{}``) Arguments to be passed to :class:`~graph_tool.inference.NestedBlockState` or :class:`~graph_tool.inference.BlockState`. bstate : :class:`~graph_tool.inference.NestedBlockState` or :class:`~graph_tool.inference.BlockState` (optional, default: ``None``) If passed, this will be used to initialize the block state directly. self_loops : bool (optional, default: ``False``) If ``True``, it is assumed that the uncertain graph can contain self-loops. References ---------- .. [peixoto-latent-2020] Tiago P. Peixoto, "Latent Poisson models for networks with heterogeneous density", Phys. Rev. E 102 012309 (2020) :doi:`10.1103/PhysRevE.102.012309`, :arxiv:`2002.07803` """def__init__(self,g,aE=numpy.nan,nested=True,state_args={},bstate=None,self_loops=False,**kwargs):super().__init__(g,nested=nested,state_args=state_args,bstate=bstate,self_loops=self_loops,**kwargs)self.q=self.g.new_ep("double",val=numpy.inf)self.q_default=-numpy.infself.S_const=0self.aE=aEifnumpy.isnan(aE):self.E_prior=Falseelse:self.E_prior=Trueself._state=libinference.make_uncertain_state(self.bstate._state,self)
[docs]defcopy(self,**kwargs):"""Return a copy of the state."""returnLatentMultigraphBlockState(**dict(self.__getstate__(),**kwargs))
def__copy__(self):returnself.copy()def__repr__(self):return"<LatentMultigraphBlockState object with %s, at 0x%x>"% \
(self.nbstateifself.nbstateisnotNoneelseself.bstate,id(self))def_mcmc_sweep(self,mcmc_state):mcmc_state.edges_only=Truereturnlibinference.mcmc_uncertain_sweep(mcmc_state,self._state,_get_rng())
[docs]classMeasuredBlockState(UncertainBaseState):r"""Inference state of a measured graph, using the stochastic block model as a prior. Parameters ---------- g : :class:`~graph_tool.Graph` Measured graph. n : :class:`~graph_tool.EdgePropertyMap` Edge property map of type ``int``, containing the total number of measurements for each edge. x : :class:`~graph_tool.EdgePropertyMap` Edge property map of type ``int``, containing the number of positive measurements for each edge. n_default : ``int`` (optional, default: ``1``) Total number of measurements for each non-edge. x_default : ``int`` (optional, default: ``0``) Total number of positive measurements for each non-edge. lp : ``float`` (optional, default: ``NaN``) Log-probability of missing edges (false negatives). If given as ``NaN``, it is assumed this is an unknown sampled from a Beta distribution, with hyperparameters given by ``fn_params`. Otherwise the values of ``fn_params`` are ignored. lq : ``float`` (optional, default: ``NaN``) Log-probability of spurious edges (false positives). If given as ``NaN``, it is assumed this is an unknown sampled from a Beta distribution, with hyperparameters given by ``fp_params`. Otherwise the values of ``fp_params`` are ignored. fn_params : ``dict`` (optional, default: ``dict(alpha=1, beta=1)``) Beta distribution hyperparameters for the probability of missing edges (false negatives). fp_params : ``dict`` (optional, default: ``dict(mu=1, nu=1)``) Beta distribution hyperparameters for the probability of spurious edges (false positives). aE : ``float`` (optional, default: ``NaN``) Expected total number of edges used in prior. If ``NaN``, a flat prior will be used instead. nested : ``boolean`` (optional, default: ``True``) If ``True``, a :class:`~graph_tool.inference.NestedBlockState` will be used, otherwise :class:`~graph_tool.inference.BlockState`. state_args : ``dict`` (optional, default: ``{}``) Arguments to be passed to :class:`~graph_tool.inference.NestedBlockState` or :class:`~graph_tool.inference.BlockState`. bstate : :class:`~graph_tool.inference.NestedBlockState` or :class:`~graph_tool.inference.BlockState` (optional, default: ``None``) If passed, this will be used to initialize the block state directly. self_loops : bool (optional, default: ``False``) If ``True``, it is assumed that the uncertain graph can contain self-loops. References ---------- .. [peixoto-reconstructing-2018] Tiago P. Peixoto, "Reconstructing networks with unknown and heterogeneous errors", Phys. Rev. X 8 041011 (2018). :doi:`10.1103/PhysRevX.8.041011`, :arxiv:`1806.07956` """def__init__(self,g,n,x,n_default=1,x_default=0,lp=numpy.nan,lq=numpy.nan,fn_params=dict(alpha=1,beta=1),fp_params=dict(mu=1,nu=1),aE=numpy.nan,nested=True,state_args={},bstate=None,self_loops=False,**kwargs):super().__init__(g,nested=nested,state_args=state_args,bstate=bstate,**kwargs)self.aE=aEifnumpy.isnan(aE):self.E_prior=Falseelse:self.E_prior=Trueself.n=nself.x=xself.n_default=n_defaultself.x_default=x_defaultself.alpha=fn_params.get("alpha",1)self.beta=fn_params.get("beta",1)self.mu=fp_params.get("mu",1)self.nu=fp_params.get("nu",1)self.lp=lpself.lq=lqself._state=libinference.make_measured_state(self.bstate._state,self)def__getstate__(self):state=super().__getstate__()returndict(state,n=self.n,x=self.x,n_default=self.n_default,x_default=self.x_default,fn_params=dict(alpha=self.alpha,beta=self.beta),fp_params=dict(mu=self.mu,nu=self.nu),lp=self.lp,lq=self.lq,aE=self.aE)
[docs]defcopy(self,**kwargs):"""Return a copy of the state."""returnMeasuredBlockState(**dict(self.__getstate__(),**kwargs))
def__repr__(self):return"<MeasuredBlockState object with %s, at 0x%x>"% \
(self.nbstateifself.nbstateisnotNoneelseself.bstate,id(self))def_mcmc_sweep(self,mcmc_state):returnlibinference.mcmc_measured_sweep(mcmc_state,self._state,_get_rng())
[docs]defset_hparams(self,alpha,beta,mu,nu):"""Set edge and non-edge hyperparameters."""self._state.set_hparams(alpha,beta,mu,nu)self.alpha=alphaself.beta=betaself.mu=muself.nu=nu
[docs]defget_p_posterior(self):"""Get beta distribution parameters for the posterior probability of missing edges."""T=self._state.get_T()M=self._state.get_M()returnM-T+self.alpha,T+self.beta
[docs]defget_q_posterior(self):"""Get beta distribution parameters for the posterior probability of spurious edges."""N=self._state.get_N()X=self._state.get_X()T=self._state.get_T()M=self._state.get_M()returnX-T+self.mu,N-X-(M-T)+self.nu
[docs]classMixedMeasuredBlockState(UncertainBaseState):r"""Inference state of a measured graph with heterogeneous errors, using the stochastic block model as a prior. Parameters ---------- g : :class:`~graph_tool.Graph` Measured graph. n : :class:`~graph_tool.EdgePropertyMap` Edge property map of type ``int``, containing the total number of measurements for each edge. x : :class:`~graph_tool.EdgePropertyMap` Edge property map of type ``int``, containing the number of positive measurements for each edge. n_default : ``int`` (optional, default: ``1``) Total number of measurements for each non-edge. x_default : ``int`` (optional, default: ``1``) Total number of positive measurements for each non-edge. fn_params : ``dict`` (optional, default: ``dict(alpha=1, beta=10)``) Beta distribution hyperparameters for the probability of missing edges (false negatives). fp_params : ``dict`` (optional, default: ``dict(mu=1, nu=10)``) Beta distribution hyperparameters for the probability of spurious edges (false positives). aE : ``float`` (optional, default: ``NaN``) Expected total number of edges used in prior. If ``NaN``, a flat prior will be used instead. nested : ``boolean`` (optional, default: ``True``) If ``True``, a :class:`~graph_tool.inference.NestedBlockState` will be used, otherwise :class:`~graph_tool.inference.BlockState`. state_args : ``dict`` (optional, default: ``{}``) Arguments to be passed to :class:`~graph_tool.inference.NestedBlockState` or :class:`~graph_tool.inference.BlockState`. bstate : :class:`~graph_tool.inference.NestedBlockState` or :class:`~graph_tool.inference.BlockState` (optional, default: ``None``) If passed, this will be used to initialize the block state directly. self_loops : bool (optional, default: ``False``) If ``True``, it is assumed that the uncertain graph can contain self-loops. References ---------- .. [peixoto-reconstructing-2018] Tiago P. Peixoto, "Reconstructing networks with unknown and heterogeneous errors", Phys. Rev. X 8 041011 (2018). :doi:`10.1103/PhysRevX.8.041011`, :arxiv:`1806.07956` """def__init__(self,g,n,x,n_default=1,x_default=0,fn_params=dict(alpha=1,beta=10),fp_params=dict(mu=1,nu=10),aE=numpy.nan,nested=True,state_args={},bstate=None,self_loops=False,**kwargs):super().__init__(g,nested=nested,state_args=state_args,bstate=bstate,**kwargs)self.aE=aEifnumpy.isnan(aE):self.E_prior=Falseelse:self.E_prior=Trueself.n=nself.x=xself.n_default=n_defaultself.x_default=x_defaultself.alpha=fn_params.get("alpha",1)self.beta=fn_params.get("beta",10)self.mu=fp_params.get("mu",1)self.nu=fp_params.get("nu",10)self._state=Noneself.q=self.g.new_ep("double")self.sync_q()self._state=libinference.make_uncertain_state(self.bstate._state,self)
[docs]defset_hparams(self,alpha,beta,mu,nu):"""Set edge and non-edge hyperparameters."""self.alpha=alphaself.beta=betaself.mu=muself.nu=nuself.sync_q()
[docs]defcopy(self,**kwargs):"""Return a copy of the state."""returnMixedMeasuredBlockState(**dict(self.__getstate__(),**kwargs))
def__copy__(self):returnself.copy()def__setstate__(self,state):self.__init__(**state)def__repr__(self):return"<MixedMeasuredBlockState object with %s, at 0x%x>"% \
(self.nbstateifself.nbstateisnotNoneelseself.bstate,id(self))
[docs]defmcmc_sweep(self,r=.5,h=.1,hstep=1,multiflip=True,**kwargs):r"""Perform sweeps of a Metropolis-Hastings acceptance-rejection sampling MCMC to sample network partitions and latent edges. The parameter ``r`` controls the probability with which edge move will be attempted, instead of partition moves. The parameter ``h`` controls the relative probability with which hyperparamters moves will be attempted, and ``hstep`` is the size of the step. The remaining keyword parameters will be passed to :meth:`~graph_tool.inference.BlockState.mcmc_sweep` or :meth:`~graph_tool.inference.BlockState.multiflip_mcmc_sweep`, if ``multiflip=True``. """returnsuper().mcmc_sweep(r=r,multiflip=multiflip,h=h,hstep=hstep,**kwargs)
[docs]classDynamicsBlockStateBase(UncertainBaseState):def__init__(self,g,s,t,x=None,aE=numpy.nan,nested=True,state_args={},bstate=None,self_loops=False,**kwargs):r"""Base state for network reconstruction based on dynamical data, using the stochastic block model as a prior. This class is not meant to be instantiated directly, only indirectly via one of its subclasses."""super().__init__(g,nested=nested,state_args=state_args,bstate=bstate,self_loops=self_loops)self.s=[g.own_property(x)forxins]self.t=[g.own_property(x)forxint]ifxisNone:x=self.u.new_ep("double")else:x=self.u.copy_property(x,g=x.get_graph())self.x=xself.aE=aEifnumpy.isnan(aE):self.E_prior=Falseelse:self.E_prior=Trueforkinkwargs.keys():v=kwargs[k]ifisinstance(v,PropertyMap):kwargs[k]=g.own_property(v)elif(isinstance(v,collections.abc.Iterable)andlen(v)>0andisinstance(v[0],PropertyMap)):kwargs[k]=[g.own_property(x)forxinv]self.params=kwargsself.os=[ns._get_any()fornsins]self.ot=[nt._get_any()forntint]self._state=self._make_state()
[docs]defset_params(self,params):r"""Sets the model parameters via the dictionary ``params``."""self.params=dict(self.params,**params)self._state.set_params(self.params)
[docs]defcopy(self,**kwargs):"""Return a copy of the state."""returntype(self)(**dict(self.__getstate__(),**kwargs))
def__copy__(self):returnself.copy()def__repr__(self):return"<%s object with %s, at 0x%x>"% \
(self.__class__.__name__,self.nbstateifself.nbstateisnotNoneelseself.bstate,id(self))def_move_proposal(self,name,beta,step,rg,transform,entropy_args):x=x_orig=self.params[name]ifisinstance(x,collections.abc.Iterable):idx=numpy.random.randint(len(x))x=x[idx]else:idx=NoneiftransformisnotNone:x=transform[1](x)ifrgisnotNone:mi=max((rg[0],x-step))ma=min((rg[1],x+step))else:mi=x-stepma=x+stepnx=numpy.random.random()*(ma-mi)+mia=0ifrgisnotNone:a-=-log(ma-mi)mi=max((rg[0],nx-step))ma=min((rg[1],nx+step))a+=-log(ma-mi)iftransformisnotNone:nx=transform[0](nx)latent_edges=entropy_args.get("latent_edges",True)density=FalseSb=self._state.entropy(latent_edges,density)ifidxisnotNone:y=self.params[name]y[idx]=nxnx=yself.set_params({name:nx})Sa=self._state.entropy(latent_edges,density)a+=beta*(Sb-Sa)ifa>0ornumpy.random.random()<exp(a):self.set_params({name:nx})returnSa-Sb,1,1ifidxisnotNone:y=self.params[name]y[idx]=x_origx_orig=yself.set_params({name:x_orig})return0,0,0
[docs]defget_edge_prob(self,u,v,x,entropy_args={},epsilon=1e-8):r"""Return conditional posterior log-probability of edge :math:`(u,v)`."""ea=get_uentropy_args(entropy_args)returnself._state.get_edge_prob(u,v,x,ea,epsilon)
[docs]defcollect_marginal(self,g=None):r"""Collect marginal inferred network during MCMC runs. Parameters ---------- g : :class:`~graph_tool.Graph` (optional, default: ``None``) Previous marginal graph. Returns ------- g : :class:`~graph_tool.Graph` New marginal graph, with internal edge :class:`~graph_tool.EdgePropertyMap` ``"eprob"``, containing the marginal probabilities for each edge. Notes ----- The posterior marginal probability of an edge :math:`(i,j)` is defined as .. math:: \pi_{ij} = \sum_{\boldsymbol A}A_{ij}P(\boldsymbol A|\boldsymbol D) where :math:`P(\boldsymbol A|\boldsymbol D)` is the posterior probability given the data. """ifgisNone:g=Graph(directed=self.g.is_directed())g.add_vertex(self.g.num_vertices())g.gp.count=g.new_gp("int",0)g.ep.count=g.new_ep("int")g.ep.eprob=g.new_ep("double")if"x"noting.ep:g.ep.xsum=g.new_ep("double")g.ep.x2sum=g.new_ep("double")g.ep.x=g.new_ep("double")g.ep.xdev=g.new_ep("double")u=self.get_graph()x=self.get_x()libinference.collect_xmarginal(g._Graph__graph,u._Graph__graph,_prop("e",u,x),_prop("e",g,g.ep.count),_prop("e",g,g.ep.xsum),_prop("e",g,g.ep.x2sum))g.gp.count+=1g.ep.eprob.fa=g.ep.count.fa/g.gp.countg.ep.x.fa=g.ep.xsum.fa/g.gp.countg.ep.xdev.fa=sqrt(g.ep.x2sum.fa/g.gp.count-g.ep.x.fa**2)returng
[docs]classEpidemicsBlockState(DynamicsBlockStateBase):def__init__(self,g,s,beta,r,r_v=None,global_beta=None,active=None,t=[],exposed=False,aE=numpy.nan,nested=True,state_args={},bstate=None,self_loops=False,**kwargs):r"""Inference state for network reconstruction based on epidemic dynamics, using the stochastic block model as a prior. Parameters ---------- g : :class:`~graph_tool.Graph` Initial graph state. s : ``list`` of :class:`~graph_tool.VertexPropertyMap` Collection of time-series with node states over time. Each entry in this list must be a :class:`~graph_tool.VertexPropertyMap` with type ``vector<int>`` containing the states of each node in each time step. A value of ``1`` means infected and ``0`` susceptible. Other values are allowed (e.g. for recovered), but their actual value is unimportant for reconstruction. If the parameter ``t`` below is given, each property map value for a given node should contain only the states for the same points in time given by that parameter. beta : ``float`` or :class:`~graph_tool.EdgePropertyMap` Initial value of the global or local transmission probability for each edge. r : ``float`` Spontaneous infection probability. r_v : :class:`~graph_tool.VertexPropertyMap` (optional, default: ``None``) If given, this will set the initial spontaneous infection probability for each node, and trigger the use of a model where this quantity is in principle different for each node. global_beta : ``float`` (optional, default: ``None``) If provided, and ``beta is None`` this will trigger the use of a model where all transmission probabilities on edges are the same, and given (initially) by this value. t : ``list`` of :class:`~graph_tool.VertexPropertyMap` (optional, default: ``[]``) If nonempty, this allows for a compressed representation of the time-series parameter ``s``, corresponding only to points in time where the state of each node changes. Each entry in this list must be a :class:`~graph_tool.VertexPropertyMap` with type ``vector<int>`` containing the points in time where the state of each node change. The corresponding state of the nodes at these times are given by parameter ``s``. active : ``list`` of :class:`~graph_tool.VertexPropertyMap` (optional, default: ``None``) If given, this specifies the points in time where each node is "active", and prepared to change its state according to the state of its neighbors. Each entry in this list must be a :class:`~graph_tool.VertexPropertyMap` with type ``vector<int>`` containing the states of each node in each time step. A value of ``1`` means active and ``0`` inactive. exposed : ``boolean`` (optional, default: ``False``) If ``True``, the data is supposed to come from a SEI, SEIR, etc. model, where a susceptible node (valued ``0``) first transits to an exposed state (valued ``-1``) upon transmission, before transiting to the infective state (valued ``1``). aE : ``float`` (optional, default: ``NaN``) Expected total number of edges used in prior. If ``NaN``, a flat prior will be used instead. nested : ``boolean`` (optional, default: ``True``) If ``True``, a :class:`~graph_tool.inference.NestedBlockState` will be used, otherwise :class:`~graph_tool.inference.BlockState`. state_args : ``dict`` (optional, default: ``{}``) Arguments to be passed to :class:`~graph_tool.inference.NestedBlockState` or :class:`~graph_tool.inference.BlockState`. bstate : :class:`~graph_tool.inference.NestedBlockState` or :class:`~graph_tool.inference.BlockState` (optional, default: ``None``) If passed, this will be used to initialize the block state directly. self_loops : bool (optional, default: ``False``) If ``True``, it is assumed that the inferred graph can contain self-loops. References ---------- .. [peixoto-network-2019] Tiago P. Peixoto, "Network reconstruction and community detection from dynamics", Phys. Rev. Lett. 123 128301 (2019), :doi:`10.1103/PhysRevLett.123.128301`, :arxiv:`1903.10833` """ifbetaisNone:beta=global_betax=kwargs.pop("x",None)ifxisNone:ifbstateisnotNone:ifnested:u=bstate.levels[0].gelse:u=bstate.gelse:u=gifisinstance(beta,PropertyMap):x=u.own_property(beta.copy())else:x=u.new_ep("double",val=beta)x.fa=log1p(-x.fa)super().__init__(g,s=s,t=t,global_beta=global_beta,active=active,x=x,r=r,r_v=r_v,exposed=exposed,aE=aE,nested=nested,state_args=state_args,bstate=bstate,self_loops=self_loops,**kwargs)def_make_state(self):returnlibinference.make_epidemics_state(self.bstate._state,self)def__setstate__(self,state):beta=state.pop("beta",None)if"x"notinstate:state["global_beta"]=betaself.__init__(**dict(state,beta=beta))
[docs]defcopy(self,**kwargs):"""Return a copy of the state."""returntype(self)(**dict(self.__getstate__(),**dict(kwargs,beta=kwargs.get("beta",None))))
[docs]defset_params(self,params):r"""Sets the model parameters via the dictionary ``params``."""self.params=dict(self.params,**params)self._state.set_params(self.params)beta=self.params["global_beta"]ifbetaisnotNone:self.x.fa=log1p(-beta)self._state.reset_m()
[docs]defmcmc_sweep(self,r=.5,p=.1,pstep=.1,h=.1,hstep=1,xstep=.1,multiflip=True,**kwargs):r"""Perform sweeps of a Metropolis-Hastings acceptance-rejection sampling MCMC to sample network partitions and latent edges. The parameter ``r`` controls the probability with which edge move will be attempted, instead of partition moves. The parameter ``h`` controls the relative probability with which moves for the parameters ``r_v`` will be attempted, and ``hstep`` is the size of the step. The parameter ``p`` controls the relative probability with which moves for the parameters ``global_beta`` and ``r`` will be attempted, and ``pstep`` is the size of the step. The paramter ``xstep`` determines the size of the attempted steps for the edge transmission probabilities. The remaining keyword parameters will be passed to :meth:`~graph_tool.inference.BlockState.mcmc_sweep` or :meth:`~graph_tool.inference.BlockState.multiflip_mcmc_sweep`, if ``multiflip=True``. """returnsuper().mcmc_sweep(r=r,p=p,pstep=p,h=h,hstep=hstep,xstep=xstep,multiflip=multiflip,**kwargs)
[docs]classIsingBaseBlockState(DynamicsBlockStateBase):def__init__(self,g,s,beta,x=None,h=None,t=None,aE=numpy.nan,nested=True,state_args={},bstate=None,self_loops=False,has_zero=False,**kwargs):r"""Base state for network reconstruction based on the Ising model, using the stochastic block model as a prior. This class is not supposed to be instantiated directly. Instead one of its specialized subclasses must be used, which have the same signature: :class:`IsingGlauberBlockState`, :class:`PseudoIsingBlockState`, :class:`CIsingGlauberBlockState`, :class:`PseudoCIsingBlockState`. Parameters ---------- g : :class:`~graph_tool.Graph` Initial graph state. s : ``list`` of :class:`~graph_tool.VertexPropertyMap` or :class:`~graph_tool.VertexPropertyMap` Collection of time-series with node states over time, or a single time-series. Each time-series must be a :class:`~graph_tool.VertexPropertyMap` with type ``vector<int>`` containing the Ising states (``-1`` or ``+1``) of each node in each time step. If the parameter ``t`` below is given, each property map value for a given node should contain only the states for the same points in time given by that parameter. beta : ``float`` Initial value of the global inverse temperature. x : :class:`~graph_tool.EdgePropertyMap` (optional, default: ``None``) Initial value of the local coupling for each edge. If not given, a uniform value of ``1`` will be used. h : :class:`~graph_tool.VertexPropertyMap` (optional, default: ``None``) If given, this will set the initial local fields of each node. Otherwise a value of ``0`` will be used. t : ``list`` of :class:`~graph_tool.VertexPropertyMap` (optional, default: ``[]``) If nonempty, this allows for a compressed representation of the time-series parameter ``s``, corresponding only to points in time where the state of each node changes. Each entry in this list must be a :class:`~graph_tool.VertexPropertyMap` with type ``vector<int>`` containing the points in time where the state of each node change. The corresponding state of the nodes at these times are given by parameter ``s``. aE : ``float`` (optional, default: ``NaN``) Expected total number of edges used in prior. If ``NaN``, a flat prior will be used instead. nested : ``boolean`` (optional, default: ``True``) If ``True``, a :class:`~graph_tool.inference.NestedBlockState` will be used, otherwise :class:`~graph_tool.inference.BlockState`. state_args : ``dict`` (optional, default: ``{}``) Arguments to be passed to :class:`~graph_tool.inference.NestedBlockState` or :class:`~graph_tool.inference.BlockState`. bstate : :class:`~graph_tool.inference.NestedBlockState` or :class:`~graph_tool.inference.BlockState` (optional, default: ``None``) If passed, this will be used to initialize the block state directly. self_loops : bool (optional, default: ``False``) If ``True``, it is assumed that the inferred graph can contain self-loops. has_zero : bool (optional, default: ``False``) If ``True``, the three-state "Ising" model with values ``{-1,0,1}`` is used. References ---------- .. [peixoto-network-2019] Tiago P. Peixoto, "Network reconstruction and community detection from dynamics", Phys. Rev. Lett. 123 128301 (2019), :doi:`10.1103/PhysRevLett.123.128301`, :arxiv:`1903.10833` """ifisinstance(s,PropertyMap):s=[s]ifisinstance(t,PropertyMap):t=[t]iftisNone:t=[]ifhisNone:h=g.new_vp("double")ifbstateisnotNone:ifnested:u=bstate.levels[0].gelse:u=bstate.gelse:u=gifx==1:x=u.new_ep("double",1)elifxisNone:x=u.new_ep("double",numpy.random.random(g.num_edges())-1/2)super().__init__(g,s=s,t=t,beta=beta,x=x,aE=aE,nested=nested,state_args=state_args,bstate=bstate,self_loops=self_loops,h=h,has_zero=has_zero,**kwargs)def__setstate__(self,state):beta=state.pop("beta",None)if"x"notinstate:state["x"]=1self.__init__(**dict(state,beta=beta))
[docs]defmcmc_sweep(self,r=.5,p=.1,pstep=.1,h=.1,hstep=1,xstep=.1,multiflip=True,**kwargs):r"""Perform sweeps of a Metropolis-Hastings acceptance-rejection sampling MCMC to sample network partitions and latent edges. The parameter ``r`` controls the probability with which edge move will be attempted, instead of partition moves. The parameter ``h`` controls the relative probability with which moves for the parameters ``r_v`` will be attempted, and ``hstep`` is the size of the step. The parameter ``p`` controls the relative probability with which moves for the parameters ``global_beta`` and ``r`` will be attempted, and ``pstep`` is the size of the step. The paramter ``xstep`` determines the size of the attempted steps for the edge coupling parameters. The remaining keyword parameters will be passed to :meth:`~graph_tool.inference.BlockState.mcmc_sweep` or :meth:`~graph_tool.inference.BlockState.multiflip_mcmc_sweep`, if ``multiflip=True``. """returnsuper().mcmc_sweep(r=r,p=p,pstep=p,h=h,hstep=hstep,xstep=xstep,multiflip=multiflip,**kwargs)
[docs]classIsingGlauberBlockState(IsingBaseBlockState):def__init__(self,*args,**kwargs):r"""State for network reconstruction based on the Glauber dynamics of the Ising model, using the stochastic block model as a prior. See documentation for :class:`IsingBaseBlockState` for details. """super().__init__(*args,**kwargs)def_make_state(self):returnlibinference.make_ising_glauber_state(self.bstate._state,self)def_mcmc_sweep(self,mcmc_state):returnlibinference.mcmc_ising_glauber_sweep(mcmc_state,self._state,_get_rng())def_mcmc_sweep_h(self,mcmc_state):returnlibinference.mcmc_ising_glauber_sweep_h(mcmc_state,self._state,_get_rng())
[docs]classCIsingGlauberBlockState(IsingBaseBlockState):def__init__(self,*args,**kwargs):r"""State for network reconstruction based on the Glauber dynamics of the continuous Ising model, using the stochastic block model as a prior. See documentation for :class:`IsingBaseBlockState` for details. Note that in this case the ``s`` parameter should contain property maps of type ``vector<double>``, with values in the range :math:`[-1,1]`. """super().__init__(*args,**kwargs)def_make_state(self):returnlibinference.make_cising_glauber_state(self.bstate._state,self)def_mcmc_sweep(self,mcmc_state):returnlibinference.mcmc_cising_glauber_sweep(mcmc_state,self._state,_get_rng())def_mcmc_sweep_h(self,mcmc_state):returnlibinference.mcmc_cising_glauber_sweep_h(mcmc_state,self._state,_get_rng())
[docs]classPseudoIsingBlockState(IsingBaseBlockState):def__init__(self,*args,**kwargs):r"""State for network reconstruction based on the equilibrium configurations of the Ising model, using the Pseudolikelihood approximation and the stochastic block model as a prior. See documentation for :class:`IsingBaseBlockState` for details. Note that in this model "time-series" should be interpreted as a set of uncorrelated samples, not a temporal sequence. """super().__init__(*args,**kwargs)def_make_state(self):returnlibinference.make_pseudo_ising_state(self.bstate._state,self)def_mcmc_sweep(self,mcmc_state):returnlibinference.mcmc_pseudo_ising_sweep(mcmc_state,self._state,_get_rng())def_mcmc_sweep_h(self,mcmc_state):returnlibinference.mcmc_pseudo_ising_sweep_h(mcmc_state,self._state,_get_rng())
[docs]classPseudoCIsingBlockState(IsingBaseBlockState):def__init__(self,*args,**kwargs):r"""State for network reconstruction based on the equilibrium configurations of the continuous Ising model, using the Pseudolikelihood approximation and the stochastic block model as a prior. See documentation for :class:`IsingBaseBlockState` for details. Note that in this model "time-series" should be interpreted as a set of uncorrelated samples, not a temporal sequence. Additionally, the ``s`` parameter should contain property maps of type ``vector<double>``, with values in the range :math:`[-1,1]`. """super().__init__(*args,**kwargs)def_make_state(self):returnlibinference.make_pseudo_cising_state(self.bstate._state,self)def_mcmc_sweep(self,mcmc_state):returnlibinference.mcmc_pseudo_cising_sweep(mcmc_state,self._state,_get_rng())def_mcmc_sweep_h(self,mcmc_state):returnlibinference.mcmc_pseudo_cising_sweep_h(mcmc_state,self._state,_get_rng())
[docs]defmarginal_multigraph_entropy(g,ecount):r"""Compute the entropy of the marginal latent multigraph distribution. Parameters ---------- g : :class:`~graph_tool.Graph` Marginal multigraph. ecount : :class:`~graph_tool.EdgePropertyMap` Vector-valued edge property map containing the counts of edge multiplicities. Returns ------- eh : :class:`~graph_tool.EdgePropertyMap` Marginal entropy of edge multiplicities. Notes ----- The mean posterior marginal multiplicity distribution of a multi-edge :math:`(i,j)` is defined as .. math:: \pi_{ij}(w) = \sum_{\boldsymbol G}\delta_{w,G_{ij}}P(\boldsymbol G|\boldsymbol D) where :math:`P(\boldsymbol G|\boldsymbol D)` is the posterior probability of a multigraph :math:`\boldsymbol G` given the data. The corresponding entropy is therefore given (in nats) by .. math:: \mathcal{S}_{ij} = -\sum_w\pi_{ij}(w)\ln \pi_{ij}(w). """eh=g.new_ep("double")libinference.marginal_count_entropy(g._Graph__graph,_prop("e",g,ecount),_prop("e",g,eh))returneh