CoCalc Public Filesmonic-sigmoidal.sagews
Author: Namig Guliyev
Views : 51
Compute Environment: Ubuntu 18.04 (Deprecated)
import sys

def continued_fraction(r):
a, b = floor(r), r - floor(r)
value = [a]
while b != 0:
a, b = floor(1/b), 1/b - floor(1/b)
value += [a]
return value

def positive_rational(n):
""" Returns the n-th term of the Calkin-Wilf sequence """
b = bin(n)[2:]
ld = '1'
c = 0
cf = []
for i in range(len(b)-1, -1, -1):
if b[i] == ld:
c += 1
else:
cf += [c]
ld = b[i]
c = 1
cf += [c]
r = cf[-1]
for k in range(len(cf)-2, -1, -1):
r = cf[k] + 1 / r
return r

def rational_number(n):
""" Returns the n-th rational number """
if n == 0:
return 0
elif n % 2 == 0:
return positive_rational(n // 2)
else:
return -rational_number(n+1)

def polynomial(n):
""" Returns the n-th monic polynomial as a list of coefficients """
if n == 1:
return [1]
else:
cf = continued_fraction(positive_rational(n))
if len(cf) == 1:
cf = [cf[0] - 2]
else:
cf = [cf[0]] + [k-1 for k in cf[1:-1]] + [cf[-1]-2]
return [rational_number(k) for k in cf] + [1]

def polynomial_function(p):
""" list of coefficients --> function """
def u(t):
value = 0
for i in range(len(p)):
value += p[i] * (t ** i)
return value
return u

beta1 = lambda t: exp(-1/t) if t > 0 else 0

beta = lambda a, b, t: beta1(b-t) / (beta1(b-t) + beta1(t-a))  # smooth transition

h = lambda d, lambda_, t: 1 - min(0.5, lambda_) / (log(t - d + 1) + 1) if t > d else 1 - min(0.5, lambda_)

coef_lin_map = lambda a, b, c, d: ((b*c - a*d) / (b - a), (d - c) / (b - a))  # coefficients of the linear map [a,b] -> [c,d]

lower_bound = lambda p: p[0] + sum([c - abs(c) for c in p[1:]]) / 2

upper_bound = lambda p: p[0] + sum([c + abs(c) for c in p[1:]]) / 2

upper_bound_for_derivative = lambda p, a, b: sum([k * abs(p[k]) * (max(abs(a), abs(b)) ** (k-1)) for k in range(1, len(p))])

def sigma(d, lambda_):
def sigma1(t):
if t >= d:
m = round(t / (2*d))
if t == (2*m+1)*d:  # preventing possible rounding error
m += 1
p = polynomial(m)
um = polynomial_function(p)
if t <= 2*m*d:
M = h(d, lambda_, (2*m+1)*d)
if len(p) == 1:  # um is constant
value = (1 + M) / 2
else:
A1, A2 = lower_bound(p), upper_bound(p)
am, bm = coef_lin_map(A1, A2, (2*M + 1)/3, (M + 2)/3)
value = am + bm * um(t / d - 2 * m + 1)
else:
K = (sigma1(2*m*d) + sigma1((2*m+1)*d)) / 2
if t - 2*m*d < d / 2:
M = h(d, lambda_, (2*m+1)*d)
if len(p) == 1:  # um is constant
v = (1 + M) / 2
delta = d / 2
else:
A1, A2 = lower_bound(p), upper_bound(p)
am, bm = coef_lin_map(A1, A2, (2*M + 1)/3, (M + 2)/3)
v = am + bm * um(t / d - 2 * m + 1)
epsilon = (1 - M) / 6
C = upper_bound_for_derivative(p, 1, 1.5)
delta = min(d / 2, (epsilon * d) / (bm * C))
value = K - (K - v) * beta(2*m*d, 2*m*d + delta, t)
else:
p = polynomial(m+1)
um = polynomial_function(p)
M = h(d, lambda_, (2*m+3)*d)
if len(p) == 1:  # um is constant
v = (1 + M) / 2
delta = d / 2
else:
A1, A2 = lower_bound(p), upper_bound(p)
am, bm = coef_lin_map(A1, A2, (2*M + 1)/3, (M + 2)/3)
v = am + bm * um(t / d - 2 * m - 1)
epsilon = (1 - M) / 6
C = upper_bound_for_derivative(p, -0.5, 0)
delta = min(d / 2, (epsilon * d) / (bm * C))
value = K - (K - v) * (1 - beta((2*m+1)*d - delta, (2*m+1)*d, t))
else:
value = (1 + h(d, lambda_, 3*d)) * (1 - beta1(d - t)) / 2
return value
return sigma1

def anti_fusc(r):
if r == 0:
return 0
cf = continued_fraction(r)
if len(cf) % 2 == 0:
cf = cf[:-1] + [cf[-1]-1] + [1]
value = 0
for i in range(len(cf)-1, -1, -1):
try:
value = value * (2 ** cf[i])
except:     # Error: Position of the polynomial is too large to calculate
return -2
if i % 2 == 0:
value += (2 ** cf[i]) - 1
return value

def anti_rational(r):
""" Returns the index of r """
if r >= 0:
return 2 * anti_fusc(r)
else:
return 2 * anti_fusc(-r) - 1

def anti_polynomial(p_):
if p_ == [1]:
return 1
del p_[-1]
p = [anti_rational(r) for r in p_]
if len(p) == 1:
p = [p[0] + 2]
else:
p = [p[0]] + [i+1 for i in p[1:-1]] + [p[-1]+2]
r = p[-1]
for k in range(len(p)-2, -1, -1):
try:
r = p[k] + 1 / r
except (MemoryError, OverflowError):     # Error: Position of the polynomial is too large to calculate
return 0
try:
return anti_fusc(r)
except (MemoryError, OverflowError):     # Error: Position of the polynomial is too large to calculate
return 0

def rational_approximating_number(x, epsilon):
for k in range(1, ceil(1 / (2 * epsilon)) + 1):
y = round(x * k) / k
if abs(x - y) <= epsilon:
return y

def rational_approximating_Taylor_polynomial(f, a, b, epsilon):
epsilon1 = epsilon / 2
n = 0
err = plot(abs(diff(f, x)), (x, a, b)).ymax() * (b - a) / 2
while (err >= epsilon1) & (n < 101):
n += 1
err = plot(abs(diff(f, x, n+1)), (x, a, b)).ymax() * ((b - a) ** (n+1)) / (factorial(n+1) * (2 ** (n+1)))
if (n == 101):
return [0, 0]
fp = taylor(f, x, (a + b) / 2, n).coefficients(sparse=False)
tp = []
for k in range(len(fp)):
tt = 0
for m in range(k, len(fp)):
tt += fp[m] * (a ** (m-k)) * binomial(m, k)
tp += [tt * ((b - a) ** k)]
epsilon2 = epsilon - err
p = []
for m in range(len(tp)-1, -1, -1):
if m == 0:
epsilon2 *= 2
r = rational_approximating_number(tp[m].n(), epsilon2 / 2)
p = [r] + p
epsilon2 -= abs(tp[m].n() - r.n())
while (p[-1] == 0) & (len(p) > 1):
del p[-1]
return p

def coefficients_of_approximation(lambda_, f, a, b, epsilon):
d = b - a
p = rational_approximating_Taylor_polynomial(f, a, b, epsilon)
if p == [0, 0]:
return [0, 0, 0, 1]
plc = p[-1]
monic = [k / plc for k in p] + [1] if plc != 0 else [1]
m = anti_polynomial(monic)
if m <= 0:                 # Error!
return [0, 0, 0, 0]
M = h(d, lambda_, (2*m+1)*d)
if m == 1:
am, bm = 1 / 2, h(d, lambda_, 3*d) / 2
else:
A1, A2 = lower_bound(monic), upper_bound(monic)
am, bm = coef_lin_map(A1, A2, (2*M + 1)/3, (M + 2)/3)
c1, c2 = plc / bm, 2 * am * plc / (bm * (1 + h(d, lambda_, 3*d)))
return [c1, c2, b - (2 * m * d), 2 * a - b]

def graph_of_sigma(d, lambda_, xmin, xmax, text_ = ""):
s = sigma(d, lambda_)
h_ = lambda t: 1 - min(0.5, lambda_) / (log(t - d + 1) + 1) if t >= d else 1 - min(0.5, lambda_)
image = plot(s, (x, xmin, xmax), color = 'blue') + plot(h_, (x, xmin, xmax), color = 'red', linestyle = "--") + plot(1, (x, xmin, xmax), color = 'green', linestyle = "--")
if text_ == "lambda":
image += text(r"$\lambda = " + str(lambda_).rstrip('0') + "$", (0.95 * xmax, 0.5 * image.ymin() + 0.5 * image.ymax()), fontsize = 15, rgbcolor = (0, 0, 0))
elif text_ == "s":
image += text(r"$s = " + str(d).rstrip('0.') + "$", (0.95 * xmax, 0.5 * image.ymin() + 0.5 * image.ymax()), fontsize = 15, rgbcolor = (0, 0, 0))
return image

#graph_of_sigma(2.0, 0.250, 0, 20).save('monic-0-20.png')
#graphics_array([[graph_of_sigma(2.0, 0.6, 0, 100, True)], [graph_of_sigma(2.0, 0.3, 0, 100, True)], [graph_of_sigma(2.0, 0.1, 0, 100, True)]]).save('monic-0-100.png', figsize = [8, 12])
#graph_of_sigma(3.0, 0.5, 0, 50).save('monic-0-50.png')
#graphics_array([[graph_of_sigma(1.0, 0.5, 0, 100, "lambda"), graph_of_sigma(1.0, 0.2, 0, 100, "lambda")], [graph_of_sigma(1.0, 0.1, 0, 100, "lambda"), graph_of_sigma(1.0, 0.05, 0, 100, "lambda")]]).save('monic-0-100_lambda.png', figsize = [12, 6])
#graphics_array([[graph_of_sigma(1.0, 0.75, 0, 100, "s"), graph_of_sigma(3.0, 0.75, 0, 100, "s")], [graph_of_sigma(6.0, 0.75, 0, 100, "s"), graph_of_sigma(9.0, 0.75, 0, 100, "s")]]).save('monic-0-100_s.png', figsize = [12, 6])

sig = sigma(2, 0.25)
tbl = []
arg = [0.4 * i for i in range(0, 50)]
for i in range(10):
tbl += [[round(arg[i], 5), round(sig(arg[i]), 5), round(arg[10+i], 5), round(sig(arg[10+i]), 5), round(arg[20+i], 5), round(sig(arg[20+i]), 5), round(arg[30+i], 5), round(sig(arg[30+i]), 5), round(arg[40+i], 5), round(sig(arg[40+i]), 5)]]
#latex(table(tbl, frame=True))

sig2 = sigma(3, 0.5)
tbl2 = []
for i in range(10):
tbl2 += [[round(i, 5), round(sig2(i), 5), round(10+i, 5), round(sig2(10+i), 5), round(20+i, 5), round(sig2(20+i), 5), round(30+i, 5), round(sig2(30+i), 5), round(40+i, 5), round(sig2(40+i), 5)]]
#latex(table(tbl2, frame=True))

def print_table_of_parameters(f, a, b, errors):
tbl = []
for epsilon in errors:
c1, c2, theta1, theta2 = coefficients_of_approximation(1/4, f, a, b, epsilon)
tbl += [[2, c1.n(25), c2.n(25), theta1.n(25), theta2, epsilon.n()]]

#f(x) = 1 + x + x^2 / 2 + x^3 / 6 + x^4 / 24 + x^5 / 120 + x^6 / 720; print_table_of_parameters(f, -1, 1, [0.95, 0.6, 0.35, 0.1, 0.04, 0.01])
#f(x) = 4 * x / (4 + x^2); print_table_of_parameters(f, -1, 1, [0.95, 0.6, 0.35, 0.1, 0.04, 0.01])
#f(x) = sin(x) - x * cos(x + 1); print_table_of_parameters(f, -1, 1, [0.95, 0.6, 0.35, 0.1, 0.04, 0.01])

def graph_with_approximators(f, a, b, errors, colors, linestyles):
image = plot(f, (x, a, b), color = 'blue', legend_label = '$f$')
for i in range(len(errors)):
p = polynomial_function(rational_approximating_Taylor_polynomial(f, a, b, errors[i])); pol = lambda t: p((t - a) / (b - a))
image += plot(pol, (x, a, b), color = colors[i], linestyle = linestyles[i], legend_label = '$\mathcal{N}\$ ($\epsilon = ' + str(errors[i]).rstrip('0') + '$)')
return image

#f(x) = 1 + x + x^2 / 2 + x^3 / 6 + x^4 / 24 + x^5 / 120 + x^6 / 720; graph_with_approximators(f, -1, 1, [0.6, 0.35, 0.1], [(1, 0, 0), (2/3, 0, 1/3), (1/3, 0, 2/3)], [":", "-.", "--"]).save('monic-f1.png')
#f(x) = 4 * x / (4 + x^2); graph_with_approximators(f, -1, 1, [0.6, 0.35, 0.1], [(1, 0, 0), (2/3, 0, 1/3), (1/3, 0, 2/3)], [":", "-.", "--"]).save('monic-f2.png')
#f(x) = sin(x) - x * cos(x + 1); graph_with_approximators(f, -1, 1, [0.6, 0.35, 0.1], [(1, 0, 0), (2/3, 0, 1/3), (1/3, 0, 2/3)], [":", "-.", "--"]).save('monic-f3.png')

@interact
def _(a_b = range_slider(-10.0, 10.0, 1, (-1.0, 1.0), label = "[a, b]: "),
lambda_ = slider(0.01, 1.0, step_size=0.001, default=0.2500, label="λ: "),
x_range = range_slider(-10.0, 100.0, 1, (0, 20.0), label = "range: ")):
d = a_b[1] - a_b[0]
s = sigma(d, lambda_)
h = lambda t: 1 - min(0.5, lambda_) / (log(t - d + 1) + 1) if t >= d else 1 - min(0.5, lambda_)
xmin, xmax = x_range[0], x_range[1]
show(plot(s, (x, xmin, xmax), color = 'blue') + plot(h, (x, xmin, xmax), color = 'red', linestyle = "--") + plot(1, (x, xmin, xmax), color = 'green', linestyle = "--")) #, aspect_ratio = 1

@interact
def __(a_b = range_slider(-10.0, 10.0, 1, (-1.0, 1.0), label = "[a, b]: "),
lambda_ = slider(0.01, 1.0, step_size=0.001, default=0.2500, label="λ: "),
t = input_box(default = 0, label = "t: ")):
d = a_b[1] - a_b[0]
s = sigma(d, lambda_)
sys.stdout.write('σ(')
sys.stdout.write(str(t.n()))
sys.stdout.write(') = ')
sys.stdout.write(str(s(t).n()))

@interact
def ___(a_b = range_slider(-10.0, 10.0, 1, (-1.0, 1.0), label = "[a, b]: "),
lambda_ = slider(0.01, 1.0, step_size=0.001, default=0.2500, label="λ: "),
epsilon = slider(0.01, 1.0, step_size=0.01, default=1.0, label="ε: "),
f = input_box(default = sin(x), label = "f: ")):
p = rational_approximating_Taylor_polynomial(f, a_b[0], a_b[1], epsilon)
if p == [0, 0]:
show(plot(f, (x, a_b[0], a_b[1]), color = 'blue', legend_label = '$f$'))
print "None of the first 100 Taylor polynomials approximates f with the given accuracy."
else:
polf = polynomial_function(p)
pol1 = lambda t: polf((t - a_b[0]) / (a_b[1] - a_b[0]))
show(plot(f, (x, a_b[0], a_b[1]), color = 'blue', legend_label = '$f$') + plot(pol1, (x, a_b[0], a_b[1]), color = 'red', linestyle = "--", legend_label = '$\mathcal{N}$'))
c1, c2, theta1, theta2 = coefficients_of_approximation(lambda_, f, a_b[0], a_b[1], epsilon)
if theta2 < 0:
print "c₁ =", c1, ", c₂ =", c2, ", θ₁ =", theta1, ", θ₂ =", theta2
elif theta2 == 0:
print "Error: Position of the polynomial is too large to calculate."