From b226f06b2283689a12203e5b10eada9a586a3adf Mon Sep 17 00:00:00 2001 From: Cliff Kerr Date: Tue, 12 Dec 2023 14:19:25 -0800 Subject: [PATCH] perplexed merge --- stisim/networks/networks.py | 55 +++++++++++- tests/devtests/devtest_orientation.py | 105 ++++++++++++++++++++++ tests/devtests/devtest_static_networks.py | 27 ++++++ 3 files changed, 186 insertions(+), 1 deletion(-) create mode 100644 tests/devtests/devtest_orientation.py create mode 100644 tests/devtests/devtest_static_networks.py diff --git a/stisim/networks/networks.py b/stisim/networks/networks.py index 17bbbd12..812f825e 100644 --- a/stisim/networks/networks.py +++ b/stisim/networks/networks.py @@ -896,4 +896,57 @@ def add_pairs(self, mother_inds, unborn_inds, dur): self.contacts.p2 = np.concatenate([self.contacts.p2, unborn_inds]) self.contacts.beta = np.concatenate([self.contacts.beta, beta]) self.contacts.dur = np.concatenate([self.contacts.dur, dur]) - return len(mother_inds) \ No newline at end of file + return + +class static(Network): + """ + A network class of static partnerships converted from a networkx graph. There's no formation of new partnerships + and initialized partnerships only end when one of the partners dies. The networkx graph can be created outside STIsim + if population size is known. Or the graph can be created by passing a networkx generator function to STIsim. + + **Examples**:: + + # Generate a networkx graph and pass to STIsim + import networkx as nx + import stisim as ss + g = nx.scale_free_graph(n=10000) + ss.static(graph=g) + + # Pass a networkx graph generator to STIsim + ss.static(graph=nx.erdos_renyi_graph, p=0.0001) + + """ + def __init__(self, graph, **kwargs): + self.graph = graph + self.kwargs = kwargs + super().__init__() + return + + def initialize(self, sim): + popsize = sim.pars['n_agents'] + if callable(self.graph): + self.graph = self.graph(n = popsize, **self.kwargs) + self.validate_pop(popsize) + super().initialize(sim) + self.get_contacts() + return + + def validate_pop(self, popsize): + n_nodes = self.graph.number_of_nodes() + if n_nodes > popsize: + errormsg = f'Please ensure the number of nodes in graph {n_nodes} is smaller than population size {popsize}.' + raise ValueError(errormsg) + + def get_contacts(self): + p1s = [] + p2s = [] + for edge in self.graph.edges(): + p1, p2 = edge + p1s.append(p1) + p2s.append(p2) + self.contacts.p1 = np.concatenate([self.contacts.p1, p1s]) + self.contacts.p2 = np.concatenate([self.contacts.p2, p2s]) + self.contacts.beta = np.concatenate([self.contacts.beta, np.ones_like(p1s)]) + return + + diff --git a/tests/devtests/devtest_orientation.py b/tests/devtests/devtest_orientation.py new file mode 100644 index 00000000..a7d1194c --- /dev/null +++ b/tests/devtests/devtest_orientation.py @@ -0,0 +1,105 @@ +""" +Experimenting with a sexual orentation class +""" + +import sciris +import stisim as ss +import numpy as np + +######################################################### +# Method using only 1 sexual network and defining sexual orientation +######################################################### +class orientation: + def __init__(self, options=None, pop_shares=None, m_partner_prob=None): + self.options = options + self.pop_shares = pop_shares + self.m_partner_prob = m_partner_prob + + +# Motivating examples: +# 1. Sexual orientation patterns from Western countries +western_prefs = orientation( + options=[0, 1, 2, 3, 4], # Abridged Kinsey scale + m_partner_prob=[0.00, 0.005, 0.020, 0.100, 1.00], # Probability of taking a male partner + pop_shares=[0.93, 0.040, 0.005, 0.005, 0.02] + # From https://en.wikipedia.org/wiki/Demographics_of_sexual_orientation +) +orientation1 = ss.State('orientation', int, + distdict=dict(dist='choice', par1=western_prefs.m_partner_prob, + par2=western_prefs.pop_shares), + eligibility='male') + +# 2. Sexual orientation with no MSM/MF mixing +exclusive_prefs = orientation( + options=['mf', 'msm'], # No mixing + m_partner_prob=[0.0, 1.0], # Probability of taking a male partner + pop_shares=[0.95, 0.05] +) + +# 3. Sexual orientation without open MSM +closeted_prefs = orientation( + options=[0, 1, 2, 3], # last category still only partners with MSM 50% of the time + m_partner_prob=[0.0, 0.005, 0.02, 0.50], # Probability of taking a male partner + pop_shares=[0.9, 0.010, 0.03, 0.06] +) + +######################################################### +# Multidimensional preferences +######################################################### +# Make a person with an age, sex, geolocation, and preferences for all these attributes in their partner +class Person: + def __init__(self, name, sex, exact_age, geo, sex_prefs=None, age_prefs=None, geo_prefs=None): + self.name = name + self.sex = sex + self.exact_age = exact_age + self.age_bins = np.array([15,30,45,60,100]) + self.age = self.age_bins[np.digitize(self.exact_age, self.age_bins)] + self.geo = geo + self.prefs = {'sex': sex_prefs, 'age': age_prefs, 'geo': geo_prefs} + self.pop_rankings = dict() + self.pop_scores = None + self.pop_order = None + + def get_preferences(self, potentials): + # Example ranking algorithm + self.pop_scores = np.ones_like(potentials) + + for this_dim, these_prefs in p.prefs.items(): + prefs = np.zeros_like(potentials) + partner_vals = np.array([getattr(q, this_dim) for q in potentials]) + for s, sprob in these_prefs.items(): + matches = ss.true(partner_vals == s) + if len(matches): + prefs[matches] = sprob / len(matches) + self.pop_rankings[this_dim] = prefs + self.pop_scores = self.pop_scores * prefs + + self.pop_order = np.array([q.name for q in potentials])[np.argsort(self.pop_scores)][::-1] + # self.pop_order = list(self.pop_order) + + +if __name__ == '__main__': + # Run tests + + p1 = Person('p1', 1, 22, 'A', sex_prefs={0:0.8, 1:0.2}, age_prefs={15:0, 30:0.99, 45:0.01, 60:0, 100:0}, geo_prefs={'A':0.6, 'B':0.3, 'C':0.1}) + p2 = Person('p2', 0, 29, 'B', sex_prefs={0:0.0, 1:1.0}, age_prefs={15:0, 30:0.5, 45:0.5, 60:0, 100:0}, geo_prefs={'A':0.2, 'B':0.7, 'C':0.1}) + p3 = Person('p3', 0, 40, 'A', sex_prefs={0:0.0, 1:1.0}, age_prefs={15:0, 30:0.4, 45:0.6, 60:0, 100:0}, geo_prefs={'A':0.9, 'B':0.1, 'C':0.0}) + p4 = Person('p4', 1, 20, 'A', sex_prefs={0:0.02, 1:0.98}, age_prefs={15:0, 30:0.9, 45:0.1, 60:0, 100:0}, geo_prefs={'A':0.5, 'B':0.5, 'C':0.0}) + + population = [p1, p2, p3, p4] + for p in population: + potentials = [q for q in population if q != p] + p.get_preferences(potentials) + +# How does each person rank each person?? +print(f'p1: {p1.pop_order}\n') +print(f'p2: {p2.pop_order}\n') + +# These are not quite right... + +# Pairs: +# [p1, p2], [p3, p4] +# [p1, p4], [na, na] + + + diff --git a/tests/devtests/devtest_static_networks.py b/tests/devtests/devtest_static_networks.py new file mode 100644 index 00000000..0fc5fa41 --- /dev/null +++ b/tests/devtests/devtest_static_networks.py @@ -0,0 +1,27 @@ +""" +Static networks from networkx +""" + +# %% Imports and settings +import stisim as ss +import matplotlib.pyplot as plt +import networkx as nx + +ppl = ss.People(10000) + +# This example runs on one static networks + the maternal network +ppl.networks = ss.Networks( + ss.static(graph=nx.erdos_renyi_graph, p=0.0001), ss.maternal() +) + +hiv = ss.HIV() +hiv.pars['beta'] = {'static': [0.0008, 0.0008], 'maternal': [0.2, 0]} + +sim = ss.Sim(people=ppl, demographics=ss.Pregnancy(), diseases=[hiv, ss.Gonorrhea()]) +sim.initialize() +sim.run() + +plt.figure() +plt.plot(sim.tivec, sim.results.hiv.n_infected) +plt.title('HIV number of infections') +plt.show()