human
Can you provide a the Python program you requested given the following:
I am working on the file "Braid Examples.ipynb".
The first task is:
Now let us examine the behavior of cohorts as we vary the target difficulty.
This can be fairly time consuming.  We're essentially going to simulate an entire network with 25 nodes for many millions of beads, in order to get enough resolution on the graph.  We use joblib to parallelize creation of the following graphs using a map/reduce algorithm.  Each job will create ncohorts/cpu_count cohorts.
The relative vertical error in this graph is $\frac{1}{\sqrt{N}}$ with $N$ cohorts as data points, so for a very smooth graph with error ~1% we need $N=10000$.  This takes a couple days on a beefy computer.  Instead let's just do $N=100$.  The resulting curve will be noisy but it's enough to extract the salient parts.  Determining the cohorts is roughly $\mathcal{O}(B_C^2)$ in the number of beads per cohort $B_C$, so we will stop adding data points when the computation time rises for large $B_C$.
Go get a coffee now...this takes a few minutes.
Can you write Python program for it?
assistant
#Here is the Python program you requested.

dctimes = [] # Initialize our results

def job(n, difficulty, ncohorts, njob):
    begin = time.process_time()
    np.random.seed((int(time.time())*njob)%(2**32-1))# We have to seed the random number generator with something
    n.reset(target=difficulty*(1<<255))              # different for each job or they will all generate exactly 
    times = []                                       # the same results!
    lastcohort = frozenset({n.nodes[0].braids[0].vertex(0)})
    while len(list(n.nodes[0].braids[0].cohorts())) < ncohorts:  # 1% error = sqrt(N)/N
        for dummy in range(100): n.tick()
        for c in n.nodes[0].braids[0].cohorts(lastcohort): pass
        lastcohort = c
    b = n.nodes[0].braids[0]
    times.append(b.cohort_time())
    #print(difficulty, times, b.ncohorts)
    return((difficulty,numpy.mean([m for (m,s) in times]), len(b.beads), b.ncohorts, time.process_time()-begin))

def parmap(f, *args):
    """ Map the function f with arguments <args>, adding an extra argument which is its job number. """
    return Parallel(n_jobs=cpu_count())(delayed(f)(*args, job) for job in range(cpu_count()))

def gettimes(n, difficulty, ncohorts):
    def reduce(result):    
        return (difficulty, numpy.mean([m[1] for m in result]), sum([m[2] for m in result]), 
                sum([m[3] for m in result]), sum([m[4] for m in result]))
    return reduce(parmap(job, n, difficulty, ncohorts/cpu_count()))

print("(target difficulty,      mean cohort time, nbeads, ncohorts, CPU time)")
for difficulty in exp(arange(-8.5, -3.375, 0.0625)):
    if(any([x[0]==difficulty for x in dctimes])): continue
    dctimes.append(gettimes(n, difficulty, 100))
    print(dctimes[-1], ",")
    if(dctimes[-1][-1] > 5*60): break # If it takes longer than 5 CPU-minutes stop
                                      # re-run this block to add another row to dctimes.