Setup¶

In [1]:
import math
from math import sqrt as mS
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.ticker import PercentFormatter as PF
import numpy as np
import pandas as pd
from scipy.integrate import quad
import seaborn as sns

Constants & Functions¶

In [2]:
def integral(c) :
    """
    Return the arc length of y = c/x.
    lower bound is c, upper bound = 1.

    Parameters:
    c - constant.
    
    Returns:
    curve_length - Curve length.
    """
    
    def integrand(x) :
        return(mS(x**4+(c**2)) / (x**2))
    
    #Integral.
    curve_length = quad(integrand, c, 1)
    return(curve_length)





def weight(c) :
    
    """
    Return the product of arc length and 
    fraction of time canoodling. 

    Parameters:
    c - constant.
    
    Returns:
    The weight of a subset of the crowd couples.
    """
    
    return(c * integral(c)[0])





def getC(sig) :
    
    """
    Return the value C corresponding to the canoodling
    fraction of time for the couples who most frequelty appear. 

    Parameters:
    sig - Significant digits.
    
    Returns:
    C - Canoodling fraction
    """
    
    run = 1
    r = (0, 10)
    xy = {}
    while run <= sig :
        
        #Significant figures
        sig10 = 10**run
        #Get the weight for each value in the range.
        for c in [i/sig10 for i in range(int(sig10*r[0]), int(sig10*r[1]))] :
            xy[c] = weight(c)
        run += 1
        
        #Find the max weight.
        W = max(xy.values())
        #Identify the corresponding max C.
        C =  list(xy.keys())[list(xy.values()).index(W)]
        
        #Update range to pinpoint max.
        r = (C - 1/sig10,  C + 1/sig10)
        
        print("The max weight is " + str(W) + " at C = " + str(C))
    
    return(C)

Fiddler¶

In [3]:
def cData() :
    data = pd.DataFrame(columns = ['x',
                                   'y',
                                   'p'])

    for x in [i/100 for i in range(0, 101)] :
        for y in [j/100 for j in range(0, 101)]:
             
            #Actual time spent canoodling is the product of these values.
            p = x*y*100
            data.loc[len(data)] = [x, y, p]
    
    return(data)





data = cData()





def heatMap(data, option, C) :
    """
    Return the heat map of canoodling 

    Parameters:
    option - Option.  1 is the base heatplot, 2 is the base heatplot and hyperbolic curve.
    
    Returns:
    Plot.
    """
    
    sns.set()
    fig = plt.figure(figsize = (11, 8))
    ax = fig.add_subplot(xlim = (0, 1),
                         ylim = (0, 1))
    
    
    
    
    #Heatmap.
    heatmap = plt.scatter(x = data['x'],
                          y = data['y'],
                          c = data['p'],
                          cmap = 'viridis_r',
                          marker = "s",
                          s = 120,
                          alpha = 0.8,
                          zorder = 2)
    
   #Title setup.
    ax.set_title('Canoodling!', fontsize = 24)

    #X-axis setup.
    ax.set_xlabel("Partner X's Preffered Canoodling %", fontsize = 22)
    ax.xaxis.set_major_formatter(PF(xmax = 1))
    #Y-axis setup.
    ax.set_ylabel("Partner Y's Preffered Canoodling %", fontsize = 22)
    ax.yaxis.set_major_formatter(PF(xmax = 1))
    ax.tick_params(axis = 'both', which = 'major', labelsize = 16)

    #Colorbar.
    cb = plt.colorbar(heatmap,
                      format = '%.0f%%')
    
    cb.set_label('Couple Time Spent Canoodling',
                 labelpad = -100,
                 rotation = 90,
                 fontsize = 22)
    cb.ax.tick_params(labelsize = 16)
    
    #Base heatmap.
    if option == 1 :
        pass
    
    #Add on hyperbola.    
    elif option != 1 :
        if option == 2 :
            c = 0.220
        
        if option == 3 :
            c = C
            
        c_x = [i/1000 for i in range(int(1000*c), 1001)]
        c_y = []

        for x in c_x:
            c_y.append(c/x)
        
        #Hyperbola.
        ax.plot(c_x,
                c_y,
                c = 'k',
                lw = 3,
                zorder = 3)
        
        if option == 2 :
            #Points and text.
            text = [("(C, 1)", (c, 1), 3, 16),
                    ("(√C, √C)", (mS(c), mS(c)), 3, 16),
                    ("(1, C)", (1, c), 3, 16),
                    ("y = C/x", (0.48, 0.8), 1, 26)]
        
        if option == 3 :
            text = [("C = " + str(100*C) + "%", (0.7, 0.5), 1, 26)]

        for t in text :
            plt.scatter(t[1][0],
                        t[1][1],
                        lw = 3,
                        c = 'k',
                        zorder = t[2])
            
            plt.annotate(text = t[0],
                         xy = t[1],
                         c = 'k',
                         size = t[3],
                         ha = 'right',
                         va = 'top',
                         zorder = 4)

    #Save file.
    fig.savefig("2025.08.08EC" + str(option) + ".png", bbox_inches = 'tight')
In [4]:
heatMap(data, 1, 0)
In [5]:
heatMap(data, 2, 0)
In [6]:
C = getC(12)
The max value is 0.3602994628025441 at C = 0.5
The max value is 0.36053476021967157 at C = 0.49
The max value is 0.36055265487137517 at C = 0.486
The max value is 0.3605528262947058 at C = 0.4864
The max value is 0.36055282876362205 at C = 0.48636
The max value is 0.36055282877310885 at C = 0.486357
The max value is 0.3605528287732632 at C = 0.4863573
The max value is 0.36055282877326517 at C = 0.48635733
The max value is 0.36055282877326533 at C = 0.486357338
The max value is 0.3605528287732654 at C = 0.4863573375
The max value is 0.3605528287732654 at C = 0.4863573375
The max value is 0.3605528287732654 at C = 0.4863573375
In [7]:
heatMap(data, 3, C)
Rohan Lewis¶

2025.08.11¶