Editing Openai/693b7cce-38f4-800d-92f7-5e56a467bfa7 (section)

===== 1. Save the script below into knowledge_abm.py (or paste into a Jupyter/Colab cell). =====
# Install dependencies if needed:

<syntaxhighlight>pip install numpy networkx matplotlib
=== optional for saving GIFs: ===
pip install pillow

</syntaxhighlight>
# Run:

<syntaxhighlight>python knowledge_abm.py

</syntaxhighlight>

In a Jupyter/Colab notebook, the plots and animation will display inline. Running locally will open matplotlib windows; to save the animation you can uncomment the writer line near the bottom of the script.

===== <syntaxhighlight lang="python"># knowledge_abm.py =====
=== Expanded ABM with multi-channel knowledge, OCEAN/MBTI/Jung/Hero's Journey, and six-degree attenuation. ===
import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import random, math
from collections import defaultdict

=== ---------------- Params ---------------- ===
N = 30
k_nearest = 4
p_rewire = 0.12
max_degree_depth = 6
dt = 1.0
T = 150
seed = 42
np.random.seed(seed)
random.seed(seed)

alpha_base = 1.0
alpha_private = 0.5
mutation_err = 0.03
mutation_cre = 0.02
decay_lambda = 0.02
kappa_surprise = 0.6
CK_theta = 0.75
CK_p = 0.85
meta_depth_cap = 4

=== ---------------- Network ---------------- ===
G = nx.watts_strogatz_graph(n=N, k=k_nearest, p=p_rewire, seed=seed)
pos = nx.spring_layout(G, seed=seed)

=== precompute distances capped at max_degree_depth ===
shortest = dict(nx.all_pairs_shortest_path_length(G))
dist_matrix = np.full((N,N), np.inf)
for i in range(N):
    for j in range(N):
        d = shortest.get(i, {}).get(j, np.inf)
        if d <= max_degree_depth:
            dist_matrix[i,j] = d

neighbors_list = [list(G.neighbors(i)) for i in range(N)]

=== ---------------- Traits ---------------- ===
O = np.clip(np.random.normal(0.5,0.15,size=N),0,1)
C = np.clip(np.random.normal(0.5,0.12,size=N),0,1)
E = np.clip(np.random.normal(0.45,0.2,size=N),0,1)
A = np.clip(np.random.normal(0.5,0.15,size=N),0,1)
N_trait = np.clip(np.random.normal(0.4,0.18,size=N),0,1)

mbti_types = []
for i in range(N):
    e_i = 'E' if E[i] > 0.5 else 'I'
    mbti_types.append(e_i + random.choice(['S','N']) + random.choice(['T','F']) + random.choice(['J','P']))

archetypes_pool = ['Hero','Mentor','Shadow','Caregiver','Trickster','Explorer']
archetypes = [random.choice(archetypes_pool) for _ in range(N)]

hj_stages_pool = ['Ordinary','Call','Refusal','Crossing','Trials','Approach','Crisis','Return']
hj_stage_idx = np.random.randint(0, len(hj_stages_pool), size=N)

=== trust matrix ===
Tmat = np.zeros((N,N))
for i in range(N):
    for j in neighbors_list[i]:
        Tmat[i,j] = np.clip(0.5 + 0.5*np.random.rand(), 0.01, 0.99)
for i in range(N):
    for j in range(N):
        if i!=j and Tmat[i,j]==0 and dist_matrix[i,j] < np.inf:
            Tmat[i,j] = max(0.02, 0.5 '' (1.0 / (1.0 + dist_matrix[i,j])) '' np.random.rand())

contact_freq = np.zeros((N,N))
for i in range(N):
    for j in neighbors_list[i]:
        contact_freq[i,j] = np.clip(np.random.rand()*1.2, 0.05, 1.2)

=== ---------------- States ---------------- ===
b_public = np.clip(np.random.beta(2,8,size=N),0,1)
b_private = np.zeros(N)
seeds = np.random.choice(N, size=max(1,N//12), replace=False)
for s in seeds:
    b_private[s] = np.clip(0.8 + 0.15*np.random.rand(), 0, 1)

e_state = np.clip(np.random.normal(0.2,0.12,size=N),0,1)
n_state = np.clip(np.random.normal(0.3,0.2,size=N),0,1)
arousal = np.clip(np.random.normal(0.1,0.08,size=N),0,1)
lambda_i = np.clip(np.random.normal(decay_lambda, 0.01, size=N), 0.005, 0.08)

=== ---------------- Functions ---------------- ===
def degree_attenuation(i,j,visibility=1.0,alpha=1.2):
    d = dist_matrix[i,j]
    if np.isinf(d): return 0.0
    if d < 1: d = 1.0
    return visibility / (d**alpha)

def sharing_probability(sender, recipient):
    base = -3.2
    betaE, betaA, betaC, betaN, betaT, betaF = 2.0, 1.5, 2.0, 3.0, 2.2, 1.0
    expr = base + betaE''E[sender] + betaA''A[sender] - betaC''C[sender] + betaT''Tmat[sender,recipient]
    expr += betaN '' math.tanh(3.0 '' N_trait[sender] * arousal[sender])
    expr += betaF * min(contact_freq[sender,recipient],1.0)
    p = 1.0 / (1.0 + math.exp(-expr))
    return np.clip(p,0,1)

def mutation_sigma(sender, recipient):
    err = mutation_err '' (1 - C[sender]) '' (1 - Tmat[sender,recipient])
    cre = mutation_cre '' O[sender] '' (1 - A[sender])
    return max(1e-6, err + cre)

def truncated_message_passing(b_vec, depth=2):
    m_prev = {(i,j): b_vec[i] for i in range(N) for j in neighbors_list[i]}
    for d in range(depth):
        m_next = {}
        for i in range(N):
            for j in neighbors_list[i]:
                incoming = []
                weights = []
                for k in neighbors_list[i]:
                    if k==j: continue
                    if (k,i) in m_prev:
                        incoming.append(m_prev[(k,i)])
                        weights.append(Tmat[i,k]*contact_freq[k,i] + 1e-6)
                if incoming and sum(weights)>0:
                    avg = sum(w*m for w,m in zip(weights,incoming))/sum(weights)
                    smoothed = 0.35''b_vec[i] + 0.65''avg
                    m_next[(i,j)] = smoothed
                else:
                    m_next[(i,j)] = b_vec[i]
        m_prev = m_next
    return m_prev

def compute_ck_depth(b_vec, depth_cap=meta_depth_cap, theta=CK_theta, p_frac=CK_p):
    N_agents = len(b_vec)
    b_level = {1: b_vec.copy()}
    for n in range(2, depth_cap+1):
        msgs = truncated_message_passing(b_level[n-1], depth=2)
        b_n = np.zeros(N)
        for i in range(N):
            inc = [msgs[(k,i)] for k in neighbors_list[i] if (k,i) in msgs]
            b_n[i] = np.mean(inc) if inc else b_level[n-1][i]
        b_level[n] = b_n
    max_n = 0
    for n in range(1, depth_cap+1):
        ok = True
        for m in range(1,n+1):
            frac = np.sum(b_level[m] > theta) / N_agents
            if frac < p_frac:
                ok = False
                break
        if ok: max_n = n
        else: break
    return max_n

=== ---------------- Simulation ---------------- ===
CK_history = []
g_std_history = []
g_by_deg = defaultdict(list)
frames = []

seed_node = np.random.choice(N)
b_public[seed_node] = 0.95
for s in np.random.choice(range(N), size=max(1,N//10), replace=False):
    b_public[s] = np.clip(b_public[s] + 0.1*np.random.rand(),0,1)

for t in range(T):
    CK_history.append(compute_ck_depth(b_public, depth_cap=meta_depth_cap))
    g_std_history.append(np.std(b_public))
    for d in range(1, max_degree_depth+1):
        nodes_at_d = [i for i in range(N) if dist_matrix[seed_node,i] == d]
        g_by_deg[d].append(np.mean(b_public[nodes_at_d]) if nodes_at_d else np.nan)
    frames.append(b_public.copy())

    incoming_inf = np.zeros(N)
    incoming_w = np.zeros(N)
    incoming_priv = np.zeros(N)
    incoming_priv_w = np.zeros(N)
    incoming_e = np.zeros(N)
    incoming_e_w = np.zeros(N)
    incoming_n = np.zeros(N)
    incoming_n_w = np.zeros(N)

    for j in range(N):
        for i in range(N):
            if i==j: continue
            Sij = degree_attenuation(i,j,visibility=alpha_base,alpha=1.2)
            if Sij==0: continue
            p = sharing_probability(j,i)
            if dist_matrix[j,i] > 1:
                p = p * (1.0/(1.0 + dist_matrix[j,i]-1.0))
            if np.random.rand() < p:
                sigma = mutation_sigma(j,i)
                mutated = b_public[j] + np.random.normal(0, math.sqrt(sigma))
                mutated = np.clip(mutated, 0, 1)
                w = Sij '' Tmat[i,j] '' (contact_freq[j,i] if contact_freq[j,i]>0 else 0.1)
                incoming_inf[i] += w * mutated
                incoming_w[i] += w
                if b_private[j] > 0 and np.random.rand() < (0.15 + 0.6*A[j]):
                    pw = alpha_private '' Sij '' Tmat[i,j] * (1 - C[j])
                    incoming_priv[i] += pw * b_private[j]
                    incoming_priv_w[i] += pw
                ew = Sij '' Tmat[i,j] '' (0.5 + 0.5*E[j])
                incoming_e[i] += ew * e_state[j]
                incoming_e_w[i] += ew
                nw = Sij '' Tmat[i,j] '' (0.3 + 0.7*(1 - abs(hj_stage_idx[j]-hj_stage_idx[i])/len(hj_stages_pool)))
                incoming_n[i] += nw * n_state[j]
                incoming_n_w[i] += nw

    for i in range(N):
        if incoming_w[i] > 0:
            avg = incoming_inf[i] / incoming_w[i]
            strength = min(1.0, incoming_w[i] / (1.0 + incoming_w[i]))
            b_public[i] += dt '' ( -lambda_i[i]''b_public[i] + strength * (avg - b_public[i]) )
        else:
            b_public[i] += dt '' ( -lambda_i[i]''b_public[i] )

        if incoming_priv_w[i] > 0:
            avgp = incoming_priv[i] / incoming_priv_w[i]
            stick = 0.2 + 0.6*C[i]
            b_private[i] = np.clip((1-stick)''b_private[i] + stick''avgp, 0, 1)
        else:
            b_private[i] *= (1 - 0.004)

        if incoming_e_w[i] > 0:
            avg_e = incoming_e[i] / incoming_e_w[i]
            e_state[i] += dt '' ( -0.08 '' e_state[i] + 0.6 * (avg_e - e_state[i]) )
        else:
            e_state[i] += dt '' ( -0.08 '' e_state[i] )

        if incoming_n_w[i] > 0:
            avg_n = incoming_n[i] / incoming_n_w[i]
            n_state[i] += dt '' ( -0.03 '' n_state[i] + 0.4 * (avg_n - n_state[i]) )
        else:
            n_state[i] += dt '' ( -0.02 '' n_state[i] )

        surprise = abs((incoming_inf[i]/incoming_w[i]) - b_public[i]) if incoming_w[i]>0 else 0.0
        arousal[i] += dt '' ( -0.09''arousal[i] + kappa_surprise * surprise )
        arousal[i] = np.clip(arousal[i], 0, 1)

        b_public[i] = np.clip(b_public[i], 0, 1)
        e_state[i] = np.clip(e_state[i], 0, 1)
        n_state[i] = np.clip(n_state[i], 0, 1)

    # rare rewiring
    if np.random.rand() < 0.015:
        node = np.random.choice(N)
        neighs = list(G.neighbors(node))
        if neighs:
            old = np.random.choice(neighs)
            candidates = [x for x in range(N) if x!=node and not G.has_edge(node,x)]
            if candidates:
                new = np.random.choice(candidates)
                G.remove_edge(node, old)
                G.add_edge(node, new)
                neighbors_list = [list(G.neighbors(i)) for i in range(N)]
                shortest = dict(nx.all_pairs_shortest_path_length(G))
                for i in range(N):
                    for j in range(N):
                        d = shortest.get(i, {}).get(j, np.inf)
                        if d <= max_degree_depth:
                            dist_matrix[i,j] = d
                        else:
                            dist_matrix[i,j] = np.inf

=== ---------------- Visualization ---------------- ===
plt.figure(figsize=(8,3))
plt.plot(CK_history)
plt.xlabel("Timestep")
plt.ylabel("CK depth (bounded)")
plt.title("CK depth over time")
plt.show()

plt.figure(figsize=(8,3))
plt.plot(g_std_history, label="std(b_public)")
for d in range(1, max_degree_depth+1):
    arr = np.array(g_by_deg[d])
    if np.all(np.isnan(arr)): continue
    plt.plot(np.arange(len(arr)), arr, label=f"deg {d}")
plt.legend(loc='upper right')
plt.xlabel("Timestep")
plt.ylabel("Value")
plt.title("Diffusion stats")
plt.show()

fig, ax = plt.subplots(figsize=(6,6))
def update(frame):
    ax.clear()
    vals = frame
    cmap = plt.get_cmap('viridis')
    nx.draw_networkx_edges(G, pos, alpha=0.25, ax=ax)
    nx.draw_networkx_nodes(G, pos, node_color=vals, cmap=cmap, vmin=0, vmax=1, ax=ax)
    ax.set_title("Public belief (node color = b_public)")
    ax.axis('off')

ani = FuncAnimation(fig, update, frames=frames, interval=100, blit=False, repeat=False)
plt.show()

print("Simulation finished. Final CK depth:", CK_history[-1], "Final std(b_public):", g_std_history[-1])
for i in range(min(8,N)):
    print(f"A{i}: OCEAN={O[i]:.2f},{C[i]:.2f},{E[i]:.2f},{A[i]:.2f},{N_trait[i]:.2f} MBTI={mbti_types[i]} Archetype={archetypes[i]} HJ={hj_stages_pool[hj_stage_idx[i]]} b_pub={b_public[i]:.2f} b_priv={b_private[i]:.2f} e={e_state[i]:.2f} ar={arousal[i]:.2f}")

</syntaxhighlight>