https://www.youtube.com/watch?v=lh1qbp7LpFU

lh1qbp7LpFU.py

import networkx as nx
import matplotlib.pyplot as plt
from sympy import factorint


def f(n):
    result = factorint(n)
    result = result.items()
    result = [n*e//p for p, e in result]
    result = sum(result)
    result = result if n > 1 else 0
    return result 


g = nx.DiGraph()

results = []
for i in range(1, 1000001):
    g.add_edge(i, f(i))
    leafs = [i for i in g.nodes() if not list(g.successors(i))]
    percent_leaf = len(leafs) / len(g.nodes()) * 100
    results.append(percent_leaf)
    print(i, percent_leaf)

plt.plot(results)
plt.show()
##C = nx.condensation(g)
##nx.draw(C, node_size=1)
##plt.savefig('condensation_arithmetic_derivatives.pdf')
##plt.close()