adaptation_specialization_material/plots/capital_measures.py

import numpy as np
import pandas as pd
from scipy.stats import entropy

from sklearn.linear_model import LinearRegression

from matplotlib import pyplot as plt
from matplotlib.gridspec import GridSpec
import matplotlib
matplotlib.use("pgf")
matplotlib.rcParams.update(
    {
        "pgf.texsystem": "xelatex",
        "font.family": "serif",
        "font.serif": "Times New Roman",
        "text.usetex": True,
        "pgf.rcfonts": False,
    }
)
plt.rcParams["text.latex.preamble"].join([
        r"\usepackage{amsmath}",              
        r"\setmainfont{amssymb}",
])

import seaborn as sns

import argparse
from os.path import join as opj
import pickle

from cmdstanpy import CmdStanModel

parser = argparse.ArgumentParser()
parser.add_argument("--input")
args = parser.parse_args()

topics = pd.read_csv(opj(args.input, "topics.csv"))
junk = topics["label"].str.contains("Junk")
topics = topics[~junk]["label"].tolist()

fig, ax = plt.subplots()

n_topics = len(pd.read_csv(opj(args.input, "topics.csv")))
df = pd.read_csv(opj(args.input, "aggregate.csv"))
age = pd.read_csv(opj(args.input, "age.csv"))
stability = pd.read_csv(opj(args.input, "institutional_stability.csv"))

resources = pd.read_parquet(opj(args.input, "pooled_resources.parquet"))
df = df.merge(resources, left_on="bai", right_on="bai")

NR = np.stack(df[[f"start_{k+1}" for k in range(n_topics)]].values).astype(int)
NC = np.stack(df[[f"end_{k+1}" for k in range(n_topics)]].values).astype(int)
expertise = np.stack(df[[f"expertise_{k+1}" for k in range(n_topics)]].values)

NR = NR[:,~junk]
NC = NC[:,~junk]
S = np.stack(df["pooled_resources"])
S = S[:,~junk]
expertise = expertise[:,~junk]

x = NR/NR.sum(axis=1)[:,np.newaxis]
y = NC/NC.sum(axis=1)[:,np.newaxis]
S = S/S.sum(axis=1)[:,np.newaxis]

# R = np.array([
#     [((expertise[:,i]>expertise[:,i].mean())&(expertise[:,j]>expertise[:,j].mean())).mean()/((expertise[:,i]>expertise[:,i].mean())|(expertise[:,j]>expertise[:,j].mean())).mean() for j in range(len(topics))]
#     for i in range(len(topics))
# ])

nu = np.load(opj(args.input, "nu_expertise_symmetric.npy"))

print(nu)

df["research_diversity"] = np.exp(entropy(x, axis=1))
df["social_diversity"] = np.exp(entropy(np.stack(df["pooled_resources"]),axis=1))
df["intellectual_diversity"] = np.exp(entropy(expertise,axis=1))
# df["social_magnitude"] = np.log(1+np.stack(df["pooled_resources"]).sum(axis=1))
df["social_magnitude"] = np.stack(df["pooled_resources"]).sum(axis=1)


expertise_matrix = np.einsum("ki,kj->kij", expertise, expertise)
social_expertise_matrix = np.einsum("ki,kj->kij", S, S)
df["intellectual_stirling"] = 1-np.einsum("ij,kij->k", nu, expertise_matrix)
df["social_stirling"] = 1-np.einsum("ij,kij->k", nu, social_expertise_matrix)

df.fillna({
    "social_stirling": 0,
    "social_diversity": 0,
    "intellectual_diversity": 0
}, inplace=True)

df["excess_social_diversity"] = df["social_diversity"]-LinearRegression().fit(df[["intellectual_diversity"]], df["social_diversity"]).predict(df[["intellectual_diversity"]])
df["excess_social_stirling"] = df["social_stirling"]-LinearRegression().fit(df[["intellectual_stirling"]], df["social_stirling"]).predict(df[["intellectual_stirling"]])

brokerage = pd.read_csv(opj(args.input, "brokerage_2000_2009.csv"))
df = df.merge(brokerage, left_on="bai", right_on="bai")
print(df)
# df["brokerage"] = np.log(1+df["brokerage"])

df.fillna(0, inplace=True)

measures = ["intellectual_diversity", "intellectual_stirling", "excess_social_diversity", "excess_social_stirling", "social_magnitude",  "brokerage"]
labels = ["\\textbf{Intellectual diversity}", "Intellectual diversity (Stirling)", "\\textbf{Excess social diversity}", "Excess social diversity (Stirling)", "\\textbf{Power} (magnitude)",  "Power (brokerage)"]

R = np.zeros((len(measures), len(measures)))

for i, a in enumerate(measures):
    for j, b in enumerate(measures):
        if i == j:
            R[i,j] = 1
        else:
            R[i,j] = np.corrcoef(df[a], df[b])[0, 1]
            print(R[i,j])

fig, ax = plt.subplots(figsize=(4,3.2)) 

R[np.tril_indices(R.shape[0],k=-1)] = np.nan

sns.heatmap(
    R,
    cmap="Reds",
    vmin=0,
    vmax=1,
    xticklabels=labels,
    yticklabels=labels,
    ax=ax,
    annot=R,
    fmt=".2f",
    annot_kws={"fontsize": 6},
    square=True
)

# ax.xaxis.set_tick_params(rotation=45)
ax.yaxis.set_tick_params(rotation=0)
ax.set_xticklabels(ax.get_xticklabels(), rotation = 45, ha="right")

fig.savefig(
    opj(args.input, f"capital_measures.eps"),
    bbox_inches="tight",
)

df[["research_diversity", "intellectual_diversity", "social_diversity"]].agg(['mean', 'std']).to_csv(opj(args.input, "capital_measures.csv"))

fig = plt.figure(figsize=[6.4, 3.2])
df = df.merge(age, how="inner", left_on="bai", right_on="bai")
df = df.merge(stability, how="inner", left_on="bai", right_on="bai")

gs = GridSpec(1, 2, width_ratios=[2.25, 1])
ax_hist = fig.add_subplot(gs[0])
ax_violin = fig.add_subplot(gs[1])

df["age"] += 4

ax_hist.hist(df["age"], bins=np.arange(0, 60, 2), histtype="step")
ax_hist.set_xlabel("Academic age in 2019")
ax_hist.set_ylabel("Number of physicists")

n_stable = len(df[df["stable"]==True])
n_unstable = len(df[df["stable"]==False])

x_stable = n_stable/len(df)
x_unstable = n_unstable/len(df)

print(x_stable, x_unstable)

df["stable"] = df["stable"].map({
    False: f"Unstable $(n={n_unstable})$",
    True: f"Stable $(n={n_stable})$"
})

g = sns.violinplot(
    data=df,
    y="age",
    hue="stable",
    split=True,
    gap=.1,
    inner="quart",
    ax=ax_violin,
    alpha=0.5,
    common_norm=True,
    bw_adjust=0.75
)
# check axes and find which is have legend
sns.move_legend(ax_violin, title="Affiliation", loc='lower right', bbox_to_anchor=(1.1, 1), frameon=False)
ax_violin.set_ylabel("Academic age in 2019")

plt.subplots_adjust(wspace=0.25, hspace=0)
fig.savefig(opj(args.input, "age.pdf"), bbox_inches="tight")