Commit 9de341bb authored by Peter steiglechner's avatar Peter steiglechner
Browse files

stochastic Model and WS model

parent d7366093
# Opinion formation model with BC and noisy expression of opinions
# 09 May 2022
# Peter Steiglechner
# peter.steiglechner@gmail.com
# IMPORTS
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
import scipy.stats as stats
import time
import xarray as xr
import sys
import pandas as pd
import warnings
warnings.filterwarnings('ignore')
class ModelT:
def __init__(self, params, track_times):
"""
Create model by initialising an array of agents with opinions
and defining how these are updated in an update step.
"""
self.params = params
seed = params["seed"]
np.random.seed(seed)
# Define Model parameters
self.n_agents = params["n_agents"] # Number of agents
self.mu = params["mu"] # convergence rate
self.nu = params["nu"] # noise level
self.bc = params["bc"] # bounded confidence radius (bias)
self.initial_dist = params["initial_dist"] # Either "uniform" or a list of floats for a given distribution (e.g. six Americas) or, if track_times[0]!=0 the opinions of all agents.
self.T = params["T"]
# Create Agents
self.schedule = []
self.ids = np.arange(self.n_agents,dtype=int)
if track_times[0]==0 :
# initialise agents at t=0
if self.initial_dist== "uniform":
self.x = stats.uniform(0,1).rvs(size=self.n_agents)
else:
segments = np.linspace(0, 1, num=len(self.initial_dist)+1)[:len(self.initial_dist)]
self.x = np.random.choice(segments, size=self.n_agents, replace=True, p=self.initial_dist) + 1/len(segments) * np.random.random(size=self.n_agents)
else:
# adopt opinions given in self.initial_dist
if not len(self.initial_dist)==self.n_agents:
print("Error: track_times does not start with 0 but the initial_dist is not equal to the opinions of n agents")
return 0
self.x = self.initial_dist
self.time = track_times[0]
self.track_times = list(track_times)
self.running=1
self.all_x = np.empty([len(self.track_times), self.n_agents])#.astype("f4")
self.all_x[0,:] = self.x.astype("f4")
return
def run(self, t_simulation, verbose=False):
"""
update all agents for t_simulation times, store values if in track_times
"""
#if verbose: starttime = time.time()
lower = 0.0 # bound of opinion spcae
upper = 1.0 # bound of opinion space
# create nr of random updates upfront
randlist = np.random.random(size=(2*t_simulation+2))
# create random noises from Gaussian dist upfront (only for unbounded model)
#randlist_addnoise = np.random.normal(loc=0.0, scale=self.nu, size=(2*t_simulation+2))
# Create random update order upfront, correct for those agent pairs with the same index
orderlist = np.random.choice(self.ids, size=(int(t_simulation)+1, 2), replace=True)
n=0
while True: # remove all monologues (ag1 talking to ag1)
n+=1
monologs = np.where((orderlist[:,0] == orderlist[:,1]))
if len(monologs[0])==0:
break
if n==100:
print("problem")
break
orderlist[monologs, :] = np.random.choice(self.ids, size=(len(monologs[0]), 2))
# run through simulation
for t in range(1,t_simulation+1):
# 1. select agents
(ag1, ag2) = orderlist[t,:]
(x1,x2) = (self.x[ag1], self.x[ag2])
# 2. Create messages (and check if they are within the bounds)
while True:
testm1 = np.random.normal(loc=x1, scale=self.nu)
if testm1>=lower and testm1<=upper:
m1=testm1
break
while True:
testm2 = np.random.normal(loc=x2, scale=self.nu)
if testm2 >= lower and testm2 <= upper:
m2=testm2
break
# 3. Accept/reject message and update if accepted
p_m_ag1 = 1/(1+np.exp(1/self.T*(abs(m2 - x1) - self.bc)))
p_m_ag2 = 1/(1+np.exp(1/self.T*(abs(m1 - x2) - self.bc)))
# 4. Update if accepted
if randlist[2*t] < p_m_ag1:
self.x[ag1] = x1 + self.mu * (m2 - x1)
if randlist[2*t+1] < p_m_ag2:
self.x[ag2] = x2 + self.mu * (m1 - x2)
# 4. Store at the time steps given in track_times.
if t in self.track_times:
ind = self.track_times.index(t)
self.all_x[ind, :] = self.x.astype("f4")
#if verbose: print(int(time.time()-starttime), "sec")
return self.running
def mainT(expname, n, icname, mu_arr, seeds, track_times, T_arr, resolution="high", verbose=False):
'''
run a model simulation for (various) n, mu, seeds and for given IC and track_times and a resoultion of bc and nu values
'''
if resolution=="high":
eps_arr = np.arange(0.025,0.41, step=0.025)
eps_arr = np.insert(eps_arr, 0, 1e-3)
nu_arr = np.arange(0,0.301, 0.01)
nu_arr = np.delete(nu_arr, 0)
nu_arr = np.insert(nu_arr, 0, 5e-3)
nu_arr = np.insert(nu_arr, 0, 2e-3)
nu_arr = np.insert(nu_arr, 0, 1e-3)
nu_arr = np.insert(nu_arr, 0, 1e-10)
elif resolution=="low":
eps_arr = np.arange(0.05,0.41, step=0.05)
eps_arr = np.insert(eps_arr, 0, 1e-3)
nu_arr = np.arange(0,0.301, 0.03)
nu_arr = np.delete(nu_arr, 0)
nu_arr = np.insert(nu_arr, 0, 1e-3)
nu_arr = np.insert(nu_arr, 0, 1e-10)
sixAm = np.array([0.07,0.11,0.12,0.18,0.33,0.19])
icdict = {"uniform": "uniform", "sixAmerica": sixAm, "sixAmericaExtreme2":sixAm**2/np.sum(sixAm**2)}
params = dict(
n_agents = n,
seed = None,
mu = None,
nu = None,
bc = None,
initial_dist = icdict[icname] if type(icname)==str else icname, # stats.uniform(0,1)
T = None,
)
t_simulation = track_times[-1] - track_times[0]
if type(seeds)== int:
seeds = np.array([seeds])
result = np.empty([len(T_arr), len(mu_arr), len(eps_arr), len(nu_arr), len(seeds), len(track_times), params["n_agents"]])
if verbose: print(result.shape)
for nT, T in enumerate(T_arr):
params["T"] = T
for n1, mu in enumerate(mu_arr):
params["mu"] = mu
for n2, eps in enumerate(eps_arr):
params["bc"] = eps
for n3, nu in enumerate(nu_arr):
params["nu"] = nu
for n4, s in enumerate(seeds):
params["seed"] = int(s)
m = ModelT(params, track_times)
running = m.run(t_simulation, verbose=False)
result[nT, n1,n2,n3,n4, :, :] = m.all_x
d = xr.DataArray(
data= result,
name="x",
dims=["T", "mu", "bc", "nu", "seed", "t", "id"],
coords={
"T":T_arr,
"mu":mu_arr,
"bc":eps_arr,
"nu":nu_arr,
"seed": seeds,
"t": track_times,
"id": np.arange(params["n_agents"])
}
)
d = d.astype("f4")
d.attrs["initial_dist"] = params["initial_dist"] if "sixAmerica" in icname or icname=="uniform" else "continue"
# STORE
if len(seeds)==1:
fname = "data/model_"+expname+"{}Initial_seeds{}.ncdf".format(icname, seeds[0])
else:
fname = "data/model_"+expname+"{}Initial_seeds{}-{}.ncdf".format(icname, seeds[0], seeds[-1])
print("saved: ", fname)
d.to_netcdf(fname)
if verbose: print("finished")
return
if __name__=="__main__":
s0 = time.time()
print(sys.argv)
track_times = [0, 10000]
mu_arr = [0.5]
n = 10
seeds = np.arange(int(sys.argv[1]), int(sys.argv[2]))
resolution = "low"
icname = "uniform"
#old_data = xr.open_dataset("model-v8_data/model-v8_long_n100_mu0.5_T_05-001_uniformInitial_seeds0-99.ncdf", engine="netcdf4")
expname = resolution+"Res_time_"
T_arr = [0.1]
expname = expname + "T{}_".format(T_arr[0])
main(expname, n, icname, mu_arr, seeds, track_times, T_arr, resolution=resolution, verbose=False)
s1 = time.time()
print("{} min {} sec".format(int((s1-s0)/60), int(s1-s0)%60 ))
# Opinion formation model with BC and noisy expression of opinions
# 09 May 2022
# Peter Steiglechner
# peter.steiglechner@gmail.com
# IMPORTS
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
import scipy.stats as stats
import time
import xarray as xr
import sys
import pandas as pd
import warnings
warnings.filterwarnings('ignore')
class ModelWS:
def __init__(self, params, track_times):
"""
Create model by initialising an array of agents with opinions
and defining how these are updated in an update step.
"""
self.params = params
seed = params["seed"]
np.random.seed(seed)
# Define Model parameters
self.n_agents = params["n_agents"] # Number of agents
self.mu = params["mu"] # convergence rate
self.nu = params["nu"] # noise level
self.bc = params["bc"] # bounded confidence radius (bias)
self.initial_dist = params["initial_dist"] # Either "uniform" or a list of floats for a given distribution (e.g. six Americas) or, if track_times[0]!=0 the opinions of all agents.
self.k, self.p = (self.params["WS"]["k"],self.params["WS"]["p"])
self.G = nx.watts_strogatz_graph(self.n_agents, self.k, self.p, seed=seed)
self.E = np.array(list(self.G.edges))
# Create Agents
self.schedule = []
self.ids = np.arange(self.n_agents,dtype=int)
if track_times[0]==0 :
# initialise agents at t=0
if self.initial_dist== "uniform":
self.x = stats.uniform(0,1).rvs(size=self.n_agents)
else:
segments = np.linspace(0, 1, num=len(self.initial_dist)+1)[:len(self.initial_dist)]
self.x = np.random.choice(segments, size=self.n_agents, replace=True, p=self.initial_dist) + 1/len(segments) * np.random.random(size=self.n_agents)
else:
# adopt opinions given in self.initial_dist
if not len(self.initial_dist)==self.n_agents:
print("Error: track_times does not start with 0 but the initial_dist is not equal to the opinions of n agents")
return 0
self.x = self.initial_dist
self.time = track_times[0]
self.track_times = list(track_times)
self.running=1
self.all_x = np.empty([len(self.track_times), self.n_agents])#.astype("f4")
self.all_x[0,:] = self.x.astype("f4")
return
def run(self, t_simulation, verbose=False):
"""
update all agents for t_simulation times, store values if in track_times
"""
#if verbose: starttime = time.time()
lower = 0.0 # bound of opinion spcae
upper = 1.0 # bound of opinion space
# create nr of random updates upfront
#randlist = np.random.random(size=(2*t_simulation+2))
# create random noises from Gaussian dist upfront (only for unbounded model)
#randlist_addnoise = np.random.normal(loc=0.0, scale=self.nu, size=(2*t_simulation+2))
# Create random update order upfront, correct for those agent pairs with the same index
#orderlist = np.random.choice(self.ids, size=(int(t_simulation)+1, 2), replace=True)
orderlist = self.E[np.random.choice(np.array(range(len(self.E))), size=int(t_simulation)+1, replace=True)]
n=0
while True: # remove all monologues (ag1 talking to ag1)
n+=1
monologs = np.where((orderlist[:,0] == orderlist[:,1]))
if len(monologs[0])==0:
break
if n==100:
print("problem")
break
orderlist[monologs, :] = np.random.choice(self.ids, size=(len(monologs[0]), 2))
# run through simulation
for t in range(1,t_simulation+1):
# 1. select agents
(ag1, ag2) = orderlist[t,:]
(x1,x2) = (self.x[ag1], self.x[ag2])
# 2. Create messages (and check if they are within the bounds)
while True:
testm1 = np.random.normal(loc=x1, scale=self.nu)
if testm1>=lower and testm1<=upper:
m1=testm1
break
while True:
testm2 = np.random.normal(loc=x2, scale=self.nu)
if testm2 >= lower and testm2 <= upper:
m2=testm2
break
# 3. Accept/reject message and update if accepted
if abs(m2-x1)<=self.bc: #randlist[2*t] < p_m_ag1:
self.x[ag1] = x1 + self.mu * (m2 - x1)
if abs(m1-x2)<= self.bc: #randlist[2*t+1] < p_m_ag2:
self.x[ag2] = x2 + self.mu * (m1 - x2)
# 4. Store at the time steps given in track_times.
if t in self.track_times:
ind = self.track_times.index(t)
self.all_x[ind, :] = self.x.astype("f4")
#if verbose: print(int(time.time()-starttime), "sec")
return self.running
def mainWS(expname, n, icname, mu_arr, seeds, track_times, WSparam, resolution="high", verbose=False):
'''
run a model simulation for (various) n, mu, seeds and for given IC and track_times and a resoultion of bc and nu values
'''
if resolution=="high":
eps_arr = np.arange(0.025,0.41, step=0.025)
eps_arr = np.insert(eps_arr, 0, 1e-3)
nu_arr = np.arange(0,0.301, 0.01)
nu_arr = np.delete(nu_arr, 0)
nu_arr = np.insert(nu_arr, 0, 5e-3)
nu_arr = np.insert(nu_arr, 0, 2e-3)
nu_arr = np.insert(nu_arr, 0, 1e-3)
nu_arr = np.insert(nu_arr, 0, 1e-10)
elif resolution=="low":
eps_arr = np.arange(0.05,0.41, step=0.05)
eps_arr = np.insert(eps_arr, 0, 1e-3)
nu_arr = np.arange(0,0.301, 0.03)
nu_arr = np.delete(nu_arr, 0)
nu_arr = np.insert(nu_arr, 0, 1e-3)
nu_arr = np.insert(nu_arr, 0, 1e-10)
sixAm = np.array([0.07,0.11,0.12,0.18,0.33,0.19])
icdict = {"uniform": "uniform", "sixAmerica": sixAm, "sixAmericaExtreme2":sixAm**2/np.sum(sixAm**2)}
params = dict(
n_agents = n,
seed = None,
mu = None,
nu = None,
bc = None,
initial_dist = icdict[icname] if type(icname)==str else icname, # stats.uniform(0,1)
WS = WSparam
)
t_simulation = track_times[-1] - track_times[0]
if type(seeds)== int:
seeds = np.array([seeds])
result = np.empty([len(mu_arr), len(eps_arr), len(nu_arr), len(seeds), len(track_times), params["n_agents"]])
if verbose: print(result.shape)
for n1, mu in enumerate(mu_arr):
params["mu"] = mu
for n2, eps in enumerate(eps_arr):
params["bc"] = eps
for n3, nu in enumerate(nu_arr):
params["nu"] = nu
for n4, s in enumerate(seeds):
params["seed"] = int(s)
m = ModelWS(params, track_times)
running = m.run(t_simulation, verbose=False)
result[n1,n2,n3,n4, :, :] = m.all_x
d = xr.DataArray(
data= result,
name="x",
dims=["mu", "bc", "nu", "seed", "t", "id"],
coords={
"mu":mu_arr,
"bc":eps_arr,
"nu":nu_arr,
"seed": seeds,
"t": track_times,
"id": np.arange(params["n_agents"])
}
)
d = d.astype("f4")
d.attrs["WS-k"] = WSparam["k"]
d.attrs["WS-p"] = WSparam["p"]
d.attrs["initial_dist"] = params["initial_dist"] if "sixAmerica" in icname or icname=="uniform" else "continue"
# STORE
if len(seeds)==1:
fname = "data/model_"+expname+"{}Initial_seeds{}.ncdf".format(icname, seeds[0])
else:
fname = "data/model_"+expname+"{}Initial_seeds{}-{}.ncdf".format(icname, seeds[0], seeds[-1])
print("saved: ", fname)
d.to_netcdf(fname)
if verbose: print("finished")
return
if __name__=="__main__":
s0 = time.time()
print(sys.argv)
track_times = [0, 10000]
mu_arr = [0.5]
n = 10
seeds = np.arange(int(sys.argv[1]), int(sys.argv[2]))
resolution = "low"
icname = "uniform"
#old_data = xr.open_dataset("model-v8_data/model-v8_long_n100_mu0.5_T_05-001_uniformInitial_seeds0-99.ncdf", engine="netcdf4")
expname = resolution+"Res_time_"
WSparam = dict(k=10, p=0.1)
mainWS(expname, n, icname, mu_arr, seeds, track_times, WSparam, resolution=resolution, verbose=False)
s1 = time.time()
print("{} min {} sec".format(int((s1-s0)/60), int(s1-s0)%60 ))
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment