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
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)
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')
heatMap(data, 1, 0)
heatMap(data, 2, 0)
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
heatMap(data, 3, C)