# -*- coding: utf-8 -*-
"""
Created on Mon Jan  9 16:07:24 2017

@author: Mariusz
"""
import matplotlib.pyplot as plt
import numpy as np
import math

def myexp(x):
    try:
        y=math.exp(x)
    except OverflowError:
        y=float('inf')
    return y

def mypow(x,y):
    try:
        z=math.pow(x,y)
    except OverflowError:
        z=float('inf')
    except ValueError:
        z=float('nan')
    return z

def fpi(f, g, x):
    # fpi function
    e0 = f(x)
    i = 0
    while i < 100:
        x1 = g(x)
        if not np.isfinite(x1):
            break
        e1 = abs(f(x1))
        if e1 < 1.0e-10:
            break
        if abs(x1-x)/max(abs(x), 1.0) < 1.0e-6:
            break
        print('e1/e0 = %e' % (e1/e0))
        x  = x1
        e0 = e1
        i += 1
    return x
    
def run_fpi(f, g, dg):
    # try fpi from different starting points
    x = 0.75
    x = fpi(f, g, x)
    print('x = %e, f(x) = %e, dg(x)=%e' % (x, f(x), dg(x)))

    x = 0.1
    x = fpi(f, g, x)
    print('x = %e, f(x) = %e, dg(x)=%e' % (x, f(x), dg(x)))

    x = -1.5
    x = fpi(f, g, x)
    print('x = %e, f(x) = %e, dg(x)=%e' % (x, f(x), dg(x)))

# the function    
def f(x):
    return myexp(x - 2.) + mypow(x,3) - x
# fixed point function and its derivative
def g1(x):
    return myexp(x - 2.) + mypow(x,3)
def dg1(x):
    return myexp(x - 2.) + 3.*mypow(x,2)
# fixed point function and its derivative
def g2(x):
    return mypow(x - myexp(x - 2.), 1./3.)
def dg2(x):
    return 1./3.*mypow(x - myexp(x - 2.),-2./3.)*(1. - myexp(x - 2.))
# fixed point function and its derivative
# this one is derived from exp(x-2) + x*(x-1)*(x+1) = 0
def g3(x):
    return -1. - myexp(x - 2.)/(mypow(x,2) - x)
def dg3(x):
    return -(myexp(x - 2.)*(mypow(x,2) - x) - myexp(x - 2.)*(2.*x - 1.))/mypow(mypow(x,2) - x,2)

# plot the function to visually find the roots
X = np.linspace(-1.5, 2., 100)

plt.figure(0)
plt.plot(X, [f(z) for z in X])
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.draw()

plt.figure(1)
plt.plot(X, X, label='x')
plt.plot(X, [g1(z) for z in X], label='g1(x)') 
plt.grid(True)
plt.xlabel('x')
plt.ylabel('g1(x)')
plt.legend(loc = 'upper right')
plt.draw()

plt.figure(2)
plt.plot(X, X, label='x')
plt.plot(X, [g2(z) for z in X], label='g2(x)') 
plt.grid(True)
plt.xlabel('x')
plt.ylabel('g2(x)')
plt.legend(loc = 'upper right')
plt.draw()

plt.figure(3)
plt.plot(X, X, label='x')
plt.plot(X, [g3(z) for z in X], label='g3(x)') 
plt.grid(True)
plt.xlabel('x')
plt.ylabel('g3(x)')
plt.legend(loc = 'upper right')
plt.draw()

print('\n\nTrying g1\n')
run_fpi(f, g1, dg1)
print('\n\nTrying g2\n')
run_fpi(f, g2, dg2)
print('\n\nTrying g3\n')
run_fpi(f, g3, dg3)

plt.show()