Im trying to calculate the hill coefficient of two logistic functions, f(x) and g(x), their composition, c(x), and input it in a datatable. Additionally, I am inputting the functions themselves, their composition and the product of hill coefficients of f(x) and g(x). In order to calculate the hill coefficients I am using a numeric bisection method as well.
Here's my code:
# Imports
import numpy as np
from numpy import log as ln
from matplotlib import pyplot as plt
from random import randint
import sympy as sym
import pandas as pd
import math
# initializing data
data = {'row': [],
'f(x)': [],
'g(x)': [],
'f(g(x))': [],
'H_f': [],
'H_g': [],
'H_fg': [],
'Product of H_f and H_g': [],
'Does this prove hypothesis?': []
}
df = pd.DataFrame(data)
#True Randomization
def logrand(b):
ra = randint(1,b)
ra = float(ra/100)
ra = 10**ra
ra = round(ra)
ra = int(ra)
return ra
# Two Hill functions
num1 = 10
number = 1
for _ in range(num1):
# Params
u = 10
c1 = logrand(300)
r1 = randint(1,u)
k1 = logrand(300)
c2 = logrand(300)
r2 = randint(1,u)
k2 = logrand(300)
# function layout
funcf = '{}/({} e^(-{}x))'.format(c1, k1, r1)
funcg = '{}/({} e^(-{}x))'.format(c2, k2, r2)
funcc = '{}/({} e^(-{}({}/({} e^(-{}x)))))'.format(c1, k1, r1, c2, k2, r2)
# figure layout
plt.rcParams["figure.figsize"] = [7.50, 3.50]
plt.rcParams["figure.autolayout"] = True
# Hill Function for f(x) and g(x) and f(g(x))
def f(x):
return c1 / (k1 np.exp((-1*r1) * x))
def g(x):
return c2 / (k2 np.exp((-1*r2) * x))
def comp(x):
return c1 / (k1 np.exp((-1*r1) * (c2 / (k2 np.exp((-1*r2) * x)))))
# EC finder
def BisectionEC10(fa, a, b):
c = 1
x = np.linspace(a, b, 1000)
ystar = 0.10 * (fa(x).max() - fa(x).min())
while abs(fa(c) - ystar) > 0.000000001:
c = (a b) / 2
if fa(c) - ystar < 0:
a = c
elif fa(c) - ystar > 0:
b = c
# print('The EC10 of the function is: ',"{0:.15f}".format(c))
# print('Output of the function when evaluated at the EC10: ',fa(c))
return c
def BisectionEC90(fa, a, b):
c = 1
x = np.linspace(a, b, 1000)
ystar = 0.90 * (fa(x).max() - fa(x).min())
while abs(fa(c) - ystar) > 0.000000001:
c = (a b) / 2
if fa(c) - ystar < 0:
a = c
elif fa(c) - ystar > 0:
b = c
# print('The EC90 of the function is: ',"{0:.15f}".format(c))
# print('Output of the function when evaluated at the EC90: ',fa(c))
return c
# EC90 and EC10 for f(x), g(x) and f(g(x))
up = 20
lo = 0
# x = np.linspace[lo,up,1000]
x = 1
# x = sym.symbols('x')
EC90_1 = BisectionEC90(f, lo, up)
EC10_1 = BisectionEC10(f, lo, up)
EC90_2 = BisectionEC90(g, lo, up)
EC10_2 = BisectionEC10(g, lo, up)
EC90_3 = BisectionEC90(comp, lo, up)
EC10_3 = BisectionEC10(comp, lo, up)
# Hill Coefficient for f(x) and g(x)
H_1 = ln(81) / (ln(EC90_1 / EC10_1))
H_1 = round(H_1,4)
H_2 = ln(81) / (ln(EC90_2 / EC10_2))
H_2 = round(H_2,4)
H_3 = ln(81) / (ln(EC90_3 / EC10_3))
H_3 = round(H_3,4)
prod = float(H_1.real) * float(H_2.real)
if prod >= float(H_3.real):
answer = 'yes'
else:
answer = 'no'
prod = round(prod,4)
# adding all data to dataframe 2
data2 = {'row': [number],
'f(x)': [funcf],
'g(x)': [funcg],
'f(g(x))': [funcc],
'H_f': [H_1],
'H_g': [H_2],
'H_fg': [H_3],
'Product of H_f and H_g': [prod],
'Does this prove hypothesis?': [answer]
}
number = number 1
df = df.append(data2, ignore_index=True)
#final dataframe
print(df)
df.to_csv(r'/Users/*****/Desktop/Research/twoarctanfuncs.csv', index = False)
Now the problem that im stuck in is that the code just doesn't run, or that it has an infinite loop and I can't figure out where the problem is happening.
This is the error message when I force stop it:
Traceback (most recent call last):
File "/Users/*****/Documents/Python/NumericBisectionMethod/venv/twologisticfuncs.py", line 112, in <module>
EC90_1 = BisectionEC90(f, lo, up)
File "/Users/*****/Documents/Python/NumericBisectionMethod/venv/twologisticfuncs.py", line 95, in BisectionEC90
if fa(c) - ystar < 0:
File "/Users/*****/Documents/Python/NumericBisectionMethod/venv/twologisticfuncs.py", line 60, in f
return c1 / (k1 np.exp((-1*r1) * x))
KeyboardInterrupt
Thanks for any help!
CodePudding user response：
It looks like your binary search (bisection) isn't starting with a large enough span.
I quickly tried your code, and then added the statement
print (a, b, c)inside thewhileloop ofBisectionEC90. The output quickly converged to0 0.0 0.0and stayed there.I replaced that print statement with
print (abs(fa(c)), ystar)and the output converged to0.025157232704402517 0.00014330069262001276.
So it makes sense that the code never exits the loop while abs(fa(c) - ystar) > 0.000000001:
EDIT: it would be a good idea to add some error checking into the while loop of your Bisection functions to prevent this type of infinite loop. Maybe something like this:
if b - a < 0.000000001:
raise Exception("Failed to find solution; a={}, b={}".format(a,b))
CodePudding user response：
If you use a print statement or debugger, you will be able to detect the issue. You are stuck in a while loop, which means your condition is never satisfied. There could be multiple reasons, and to find out the exact issue, you would have to debug it. Although the most logical conclusion is that while abs(fa(c) - ystar) > 0.000000001 is not being satisfied as mentioned by @Mr. Snrub
I would also recommend using the if conditions to break the loop in such cases.