# this program adds a zero-crossing detector to try to find a 'root' insteadof a minima
# program to illustrate iteration, slightly more robust than the previous
# this version has a convergence criteria, and also adapts the step size
# note that your starting point determines which zero you find!
# We iterate to find an x, where log(x)/2 == cos(x)
# Note this program uses our own function eqn(x), that returns cos(x)-log(x)/2

# (try using a seed like '10', which will illustrate that even this simple 
#  program can have major logic BUGS)

import numpy as np
import matplotlib.pyplot as plt

# function definition for the equation you are trying to find the zero
def eqn(x):                       # You can put any function here
    return np.cos(x)-np.log(x)/2    # return cos(x) - log(x)

# you need to select a starting point (this will determine which zero you find)
seed = 10                   # starting point
seedx = seed                # keep the original seed, and vary seedx
region = [.2,14]            # region in which to look for a root
# convergence criteria, how many sig digits?
conv = 0.00000001
dx=seedx/1000               # starting step size towards (hopefully) convergence


x = np.linspace(region[0],region[1],100)   # set the 'grid' size for plotting (only)

# first plot what cos(x)-log(x) looks like (we could actually just look at the plot to see where eqq==0)
fig = plt.figure()
ax  = fig.add_subplot(111)
ax.plot(x,eqn(x),'b')       # plot residual of cos(x)-log(x), calls eqq()
ax.plot(x,np.cos(x),'g')    # plot log and cos functions
ax.plot(x,np.log(x)/2,'r')
ax.set_title('Log/2 (red) and cos of x (green), and value of cos(x)- log(x)/2 (blue)');
ax.set_xlabel('variable x')
ax.set_ylabel('values of functions')
ax.plot([0,region[1]],[0,0],':k')   # draw a zero line on plot

# now start the iteration to get the program to find the zero of eqq()
# we are trying to reduce the residual from our starting point
rold = eqn(seedx)           # error or residual from initial guess
seedx += dx                 # start heading to right
rnew = eqn(seedx)           # residual from slightly different guess

while (np.abs(dx) > np.abs(conv)):
    if (np.abs(rold) < np.abs(rnew)):  # if our guess made it worse, change direction                           
        dx = -dx/2                     # make dx oposite sign, and take smaller steps
    seedx += dx                        # most of the direction changes are right around the zero
    rold = rnew                        # we have made a step, so our old position becomes rold
    rnew = eqn(seedx)                  # we keep going in dx steps until we cross zero and the sign
	                                   # changes, we then turn around and take 1/2 size steps, and repeat


ax.plot([seedx,seedx],[-1,1],':y')  # plot a vertical line at the zero point of the function
ax.text(seedx,.25,"Local Minima")
plt.show()

# partial solution is to look for zero crossings, instead of minima
# go in both directions and look for a sign change
signf       = np.sign(eqn(seed))      # find the sign of the func at x
foundchange = 1                       # this is a flag to say we have found a solution
xp          = seed                    # we are going to look for a zero in both the positive xp
xm          = seed                    #  and in the negative direction xm
dx          = seed/1000               # step size
while (foundchange):                  # check the flag, if 0 then we are done
    if xp < region[1]:                # keep moving in the positive direction to look for a sign change
        xp += dx                      # step in the positive x direction and check the sign of eqn()
        signp = np.sign(eqn(xp))
        if signp != signf:            # if sign changes, then change flag and save result in xf
            foundchange = 0
            xf = xp
            
    if xm > region[0]:                # this goes in the negative x direction
        xm -= dx
        signp = np.sign(eqn(xm))
        if signp != signf:
            foundchange = 0 
            xf = xm

ax.plot([xf,xf],[-1,1],':k')  # plot a vertical line at the zero point of the function
ax.text(xf+dx,.25,"Zero crossing")