import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns
mpl.rcParams['figure.figsize'] = [16.0, 5.0]
plt.xkcd(length=200, randomness=1)
palette = sns.color_palette("deep")
from functools import partial
colors = {}
for i, label in enumerate(["productivity", "work", "toolsmiths", "time"]):
    colors[label] = palette[i]
import math
def sigmoid(x, a=1, b=1, c=1, d=1, e=0, f=0):
    return a / (b + c * math.exp(-d * (x + e))) + f
def display_curve(title, curve, xs=None, show_project_size=True):
    if xs is None:
        xs = list(range(100))
    plt.title(title)
    plt.plot(xs, [curve(x) for x in xs], color=colors["productivity"])
    plt.xlabel("Person-days of work")
    plt.ylabel("Productivity multiplier")
    if show_project_size:
        plt.axvline(x=100 * 100, color="black", linestyle="--", label="Project size")
        plt.legend()
    plt.show()
    
xs = [x/10 for x in range(-100, 100)]
display_curve("The default shape of sigmoid", sigmoid, xs, show_project_size=False)
def brookes_law(current_productivity, number_of_engineers):
    return current_productivity - min(.5, number_of_engineers * .01)
xs = list(range(0, 101))
plt.title("Work done per day per engineer with increasing engineers")
plt.xlabel("Number of engineers in the team")
plt.ylabel("Productivity Multiplier")
plt.plot(xs, [brookes_law(1, x) for x in xs], color=colors["productivity"])
plt.ylim(bottom=0)
plt.show()
def calculate_work(*, toolsmiths, engineers, productivity_curve, brookes_law=brookes_law, project_size=10_000, default_days=200):
    engineer_work = 0
    toolsmith_work = 0
    productivity = 1
    days = 0
    days_to_complete = None
    work_in_default_days = None
    ys = []
    
    while engineer_work < project_size or days < default_days:
        productivity = productivity_curve(toolsmith_work)
        toolsmith_work += toolsmiths * brookes_law(productivity, toolsmiths)
        engineer_work += engineers * brookes_law(productivity, engineers)
        days += 1
        
        ys.append(engineer_work)
        
        if engineer_work >= project_size and days_to_complete is None:
            days_to_complete = days
        if days == default_days:
            work_in_default_days = engineer_work
    
    return (work_in_default_days, days_to_complete, ys)
baseline_work, baseline_time, ys = calculate_work(toolsmiths=0, engineers=100, productivity_curve=lambda x: 1)
print(f"{baseline_work = }, {baseline_time = }")
def compare_work(productivity_curve, toolsmiths_counts, title):
    plt.title(title)
    
    for toolsmiths in toolsmiths_counts:
        _, _, ys = calculate_work(toolsmiths=toolsmiths, 
                                  engineers=100-toolsmiths, 
                                  productivity_curve=productivity_curve)
        plt.plot(list(range(0, 200)), ys[:200], label=f"{toolsmiths} toolsmiths")
    
    plt.xlabel("Total days")
    plt.ylabel("Work done")
    plt.ylim(bottom=0)
    plt.xlim(left=0)
    plt.legend()
    plt.show()
compare_work(
    productivity_curve=lambda x: 1 + x / 10_000,
    toolsmiths_counts=list(range(0, 100, 20)),
    title="Linear productivity curve")
           
def optimum_work(productivity_curve, brookes_law=brookes_law, total_engineers=100, project_size=10_000, default_days=200):
    max_work = None
    min_duration = None
    
    for toolsmiths in range(total_engineers):
        work, duration, _ = calculate_work(toolsmiths=toolsmiths, 
                                           engineers=total_engineers - toolsmiths, 
                                           productivity_curve=productivity_curve, 
                                           brookes_law=brookes_law,
                                           project_size=project_size,
                                           default_days=default_days)
        if max_work is None or max_work[1] < work:
            max_work = (toolsmiths, work)
        if min_duration is None or min_duration[1] > duration:
            min_duration = (toolsmiths, duration)
            
    return max_work, min_duration
max_work, min_time = optimum_work(lambda x: 1 + x / 10_000, brookes_law)
print(f"{max_work = }, {min_time = }")
xs = []
ys = []
zs = []
for toolsmiths in range(100):
    work, duration, _ = calculate_work(toolsmiths=toolsmiths, 
                                       engineers=100 - toolsmiths, 
                                       productivity_curve=lambda x: 1 + x / 10_000) 
    xs.append(toolsmiths)
    ys.append(work)
    zs.append(duration)
fig, ax1 = plt.subplots()
ax1.set_title("Exploring toolsmith allocation")
ax2 = ax1.twinx()
ax1.plot(xs, ys, label="Total work done in 200 days (person_days)", color=colors["work"])
ax1.set_ylabel("Work done (person-days)")
ax1.set_xlabel("# of toolsmiths (out of a team of 100)")
ax2.plot(xs, zs, label="Time to complete 10_000 person-days project (days)", color=colors["time"])
ax2.set_ylabel("Time to completion (days)")
ax1.hlines([max_work[1]], 0, max_work[0], linestyle="dotted", color=colors["work"])
ax1.vlines([max_work[0]], 0, max_work[1], linestyle="dotted", color=colors["work"])
ax2.hlines([min_time[1]], min_time[0], 100, linestyle="dotted", color=colors["time"])
ax2.vlines([min_time[0]], 0, min_time[1], linestyle="dotted", color=colors["time"])
ax1.set_xlim(0, 100)
ax1.set_ylim(bottom=0, top=12000)
ax2.set_ylim(bottom=0, top=600)
ax1.legend(loc="lower right")
ax2.legend(loc="lower left")
plt.show()
def explore_optimum_work(*, title, xlabel, xs, ys, zs, ws, vs):
    """
        xs: value along x-axis
        ys: work completed / baseline
        zs: toolsmith allocation for optimum
        ws: time to completion / baseline
        vs: toolsmith allocation for optimum
    """
    fig, axs = plt.subplots(2, 1)
    fig.set_size_inches(16.0, 10.0)
    
    axw = axs[0]
    axw.set_title(title)
    axw.plot(xs, ys, label="Total work done / baseline", color=colors["work"])
    axw.set_ylabel("Total work done / baseline")
    axw.set_xlabel(xlabel)
    
    axw.set_ylim(bottom=0)
    axw.legend(loc='upper left')
    
    axw2 = axw.twinx()
    axw2.plot(xs, zs, label="Toolsmiths", color=colors["toolsmiths"])
    axw2.set_ylabel("Toolsmiths")
    axw2.set_ylim(bottom=0, top=100)
    axw2.legend(loc='center right')
    
    axd = axs[1]
    axd.plot(xs, ws, label="Total time (days) / baseline", color=colors["time"])
    axd.set_ylabel("Total time (days) / baseline")
    axd.set_xlabel(xlabel)
    axd.set_ylim(bottom=0)
    axd.legend(loc='lower left')
    
    axd2 = axd.twinx()
    axd2.plot(xs, zs, label="Toolsmiths", color=colors["toolsmiths"])
    axd2.set_ylabel("Toolsmiths")
    axd2.set_ylim(bottom=0, top=100)
    axd2.legend(loc='center right')
    
    plt.tight_layout()
    plt.show()
linear_productivity = lambda x, multiplier=1: 1 + multiplier * x / 10_000
xs = []
ys = []
zs = []
ws = []
vs = []
base_work, base_days, _ = calculate_work(toolsmiths=0, engineers=100, productivity_curve=linear_productivity, brookes_law=lambda x, _: x)
for multiplier10 in range(0, 41):
    multiplier = multiplier10 / 10
    xs.append(multiplier + 1)
    work, duration = optimum_work(partial(linear_productivity, multiplier=multiplier), brookes_law=lambda x, _: x)
    ys.append(work[1] / base_work)
    zs.append(work[0])
    ws.append(duration[1] / base_days)
    vs.append(duration[0])
    
explore_optimum_work(
    title="Linear productivity curve without brooke's law", 
    xlabel="Final productivity at the end",
    xs=xs,
    ys=ys,
    zs=zs,
    ws=ws,
    vs=vs)
def exploratory_tooling_curve(x, offset=0):
    a = 49
    d = 1 / 12000
    e = -40000 + offset
    return 1 + sigmoid(x, a=a, d=d, e=e) - sigmoid(0, a=a, d=d, e=e)
    
display_curve("A complex domain for toolsmiths", exploratory_tooling_curve, list(range(0, 100_000, 10)))
base_work, base_days, _ = calculate_work(toolsmiths=0, engineers=100, productivity_curve=exploratory_tooling_curve)
xs = []
ys = []
zs = []
ws = []
vs = []
for offset in range(0, 100_000, 3000):
    xs.append(offset)
    work, duration = optimum_work(lambda x: exploratory_tooling_curve(x, offset=offset), brookes_law)
    ys.append(work[1] / base_work)
    zs.append(work[0])
    ws.append(duration[1] / base_days)
    vs.append(duration[0])
explore_optimum_work(
    title="A complex domain", 
    xlabel="Existing progress around the productivity curve",
    xs=xs,
    ys=ys,
    zs=zs,
    ws=ws,
    vs=vs)
compare_work(productivity_curve=exploratory_tooling_curve, 
             toolsmiths_counts=[0, 40, 80, 99], 
             title="Exploratory tooling curve")
def defined_tooling_curve(x, bonus=1):
    a = bonus
    d = 1 / 1000
    e = -5000
    return 1 + sigmoid(x, a=a, d=d, e=e) - sigmoid(0, a=a, d=d, e=e)
    
display_curve(f"A defined problem: max multiplier {1 + 1}", lambda x: defined_tooling_curve(x, bonus=1), list(range(0, 10_000)))
base_work, base_days, _ = calculate_work(toolsmiths=0, engineers=100, productivity_curve=defined_tooling_curve)
xs = []
ys = []
zs = []
ws = []
vs = []
for bonus10 in range(0, 100):
    bonus = bonus10/10
    xs.append(bonus + 1)
    work, duration = optimum_work(lambda x: defined_tooling_curve(x, bonus=bonus), brookes_law)
    ys.append(work[1] / base_work)
    zs.append(work[0])
    ws.append(duration[1] / base_days)
    vs.append(duration[0])
    
explore_optimum_work(
    title="A defined problem", 
    xlabel="Final productivity multiplier with 10,000 person-days of effort",
    xs=xs,
    ys=ys,
    zs=zs,
    ws=ws,
    vs=vs)
compare_work(productivity_curve=defined_tooling_curve, 
             toolsmiths_counts=[0, 40, 80, 99], 
             title="Defined tooling curve")
def tiny_tooling_curve(x, bonus=1):
    a = bonus
    d = 1 / 100
    e = -500
    return 1 + sigmoid(x, a=a, d=d, e=e) - sigmoid(0, a=a, d=d, e=e)
    
display_curve(f"A defined area: max multiplier {1 + 1}", lambda x: tiny_tooling_curve(x, bonus=1), list(range(0, 10_000)))
base_work, base_days, _ = calculate_work(toolsmiths=0, engineers=100, productivity_curve=tiny_tooling_curve)
xs = []
ys = []
zs = []
ws = []
vs = []
for bonus10 in range(0, 30):
    bonus = bonus10/10
    xs.append(bonus + 1)
    work, duration = optimum_work(lambda x: tiny_tooling_curve(x, bonus=bonus), brookes_law)
    ys.append(work[1] / base_work)
    zs.append(work[0])
    ws.append(duration[1] / base_days)
    vs.append(duration[0])
    
explore_optimum_work(
    title="A tiny problem", 
    xlabel="Final productivity multiplier with 10,000 person-days of effort",
    xs=xs,
    ys=ys,
    zs=zs,
    ws=ws,
    vs=vs)
compare_work(productivity_curve=tiny_tooling_curve, 
             toolsmiths_counts=[0, 40, 80, 99], 
             title="Tiny tooling curve")
def tooling_curve(x, size=1):
    a = 1
    d = 1 / (100 * size)
    e = -500 * size
    return 1 + sigmoid(x, a=a, d=d, e=e) - sigmoid(0, a=a, d=d, e=e)
display_curve(f"Varying the size: {5} of the curve", lambda x: tooling_curve(x, size=5), list(range(0, 10_000)))
base_work, base_days, _ = calculate_work(toolsmiths=0, engineers=100, productivity_curve=tooling_curve)
xs = []
ys = []
zs = []
ws = []
vs = []
for size in range(1, 100):
    xs.append(size / 10)
    work, duration = optimum_work(partial(tooling_curve, size=size/10), brookes_law)
    ys.append(work[1] / base_work)
    zs.append(work[0])
    ws.append(duration[1] / base_days)
    vs.append(duration[0])
    
explore_optimum_work(
    title="Differently sized problems", 
    xlabel="Size of problem compared to the actual project",
    xs=xs,
    ys=ys,
    zs=zs,
    ws=ws,
    vs=vs)
compare_work(productivity_curve=partial(tooling_curve, size=5), 
             toolsmiths_counts=range(0, 100, 20),
             title="Tooling curve with size 5")
from mpl_toolkits import mplot3d
import numpy as np
baseline_work, baseline_time, _ = calculate_work(toolsmiths=0, engineers=100, productivity_curve=lambda x: 1)
def final_productivity_curve(x, y):
    a = x - 1
    d = 1 / (1000 * y)
    e = -5000 * y
    adj = sigmoid(0, a=a, d=d, e=e)
    
    def curve(x):
        return 1 + sigmoid(x, a=a, d=d, e=e) - adj
    
    work, duration = optimum_work(curve)
    return (work[0], work[1] / baseline_work, duration[0], duration[1] / baseline_time)
samples = 20
x = np.linspace(1, 10, samples)
y = np.linspace(0.1, 3, samples)
X, Y = np.meshgrid(x, y)
Z = np.vectorize(final_productivity_curve)(X, Y)
def add_3d_plot(ax, z, title, label=None, cmap=mpl.cm.coolwarm):
    label = label or title
    ax.set_title(title)
    ax.set_xlabel("Final productivity multiplier")
    ax.set_ylabel("Size of the toolsmiths project")
    ax.set_zlabel(label)
    ax.plot_surface(X, Y, z, label=label, cmap=cmap)
    
elevation=20
    
fig = plt.figure(figsize=(20, 20))
ax1 = fig.add_subplot(2, 2, 1, projection='3d')
ax1.view_init(elev=elevation, azim=240)
add_3d_plot(ax1, Z[0], "Toolsmiths to maximize work", "Ideal number of toolsmiths", cmap=mpl.cm.copper)
ax2 = fig.add_subplot(2, 2, 2, projection='3d')
ax2.view_init(elev=elevation, azim=180)
add_3d_plot(ax2, Z[1], "Maximum possible work / baseline", cmap=mpl.cm.coolwarm_r)
ax3 = fig.add_subplot(2, 2, 3, projection='3d')
ax3.view_init(elev=elevation, azim=240)
add_3d_plot(ax3, Z[2], "Toolsmiths to minimize time", "Ideal number of toolsmiths", cmap=mpl.cm.copper)
ax4 = fig.add_subplot(2, 2, 4, projection='3d')
ax4.view_init(elev=elevation, azim=-45)
add_3d_plot(ax4, Z[3], "Minimimum possible time / baseline", cmap=mpl.cm.coolwarm)
plt.tight_layout()
plt.show()