monic-sigmoidal.sagews

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()]]
    print latex(table(tbl, header_row=True, frame=True))

#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."

Interact: please open in CoCalc

Interact: please open in CoCalc

Interact: please open in CoCalc