In another post I had discussed about
Golden Section Search method, how it works, implementing it in Python and some application examples also
presented like finding the shortest distance from a point to a line. If you don't
read the post, I strongly suggest to read it, because it will make easier to
follow this tutorial.
In this post I will discuss in more detail about the example, how to find the shortest distance from a point to a line. I will go through step by step until getting the result as shown in figure 1.
|
|
Figure 1. Golden Section Search application. Finding the shortest
distance from a point to a line |
How to Compute Distance Between Two Points
I begin this tutorial with a question, how to to find a distance between two points? A distance between two points can be calculated in many ways. The Euclidean distance is the most commonly used to calculate a straight line between two points on a Cartesian coordinate system. The Euclidean distance can be computed using the equation below. \[d=\sqrt{(x_1-x_2)^2+(y_1-y_2)^2} \]
In the Golden Section Search implementation, the distance computation will be done in each iteration for each interior point $x_1$ and $x_2$. That's why in figure 1 we can see two red lines from a determined blue point outside the line to each interior point. The iteration will continue and stop when it reaches a threshold error value. The threshold value is really small like 0.05 to make sure the difference between the two distances is really small in a long precision or even the same in a shorter precision.
Implementing The Solution in Python
Now let's do it in Python. Firstly we need to import some required libraries, such as: numpy, matplotlib, time and IPython.
#IMPORTING REQUIRED LIBRARIES %matplotlib inline import numpy as np import matplotlib.pyplot as plt from IPython.display import clear_output import time
Cause we are working with a line, we need to determine the line function. A function for a line can be expressed as: $f(x)=mx+b$. Where $x$ is input variable, $m$ is slope and $b$ is a constant.
The slope can be computed with the following equation. \[m=\frac{y_2-y_1}{x_2-x_1} \] Then $b$ can be found by substituted the value of $y_1$ and $x_1$ or $y_1$ and $x_1$ as respective $f(x)$ and $x$.
Below is the function to calculate the slope $m$ and $b$ with four input
variables $x_1$,$y_1$ and $x_2$,$y_2$ which are the start and end coordinates
of a line.
def line_func(x1,x2,y1,y2): m=(y2-y1)/(x2-x1) b=y1-(m*x1) return m,b
After getting $m$ and $b$, then we use it to compute $f(x)$ or $y$. For this we create a function to compute it.
def func_fx(m,b,x): fx=m*x+b return fx
Next we create a function to compute distance between two points using the equation above.
def distance(x1,y1,x2,y2): d=np.sqrt((x1-x2)**2+(y1-y2)**2) return d
To select a correct optimum point, we have to know the position of interior points one to another. For this we create a function which is called check_pos. If $x_1$ is greater than $x_2$, it means $x_1$ is on the right of $x_2$. For that we give a label 'right'.
def check_pos(x1,x2): if x2<x1: label='right' else: label='' return label
Then we create a function which is called update_interior. This function is used to update interior values $x_1$ and $x_2$ based on boundary point $x_u$ and $x_l$.
def update_interior(xl,xu): d=((np.sqrt(5)-1)/2)*(xu-xl) x1=xl+d x2=xu-d return x1,x2
The objective for this case is to find the minimum difference between the two
distances from a known point to each interior point $x_1$ and $x_2$. Therefore
we use find_min function. In the function we compute each distance as
in line 4 and 5 in the following code. From line 6 we do the distances
comparison. Based on the rule criteria then the boundary and interior
values are updated.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
def find_min(xl,xu,x1,x2,label): y1=func_fx(m,b,x1) y2=func_fx(m,b,x2) fx1=distance(x_point,y_point,x1,y1) fx2=distance(x_point,y_point,x2,y2) if fx2>fx1 and label=='right': xl=x2 xu=xu new_x=update_interior(xl,xu) x1=new_x[0] x2=new_x[1] xopt=x1 else: xl=xl xu=x1 new_x=update_interior(xl,xu) x1=new_x[0] x2=new_x[1] xopt=x2 return xl,xu,xopt,fx1,fx2 |
Lastly we create the golden search function with four input variables namely:
upper boundary($x_u$), lower boundary($x_l$), mode (which is 'min' because we
only searching for minimum value) and error threshold($et$). In the function
there is also plotting command at line 28-34 to produce a graph like figure 1.
The rest of this function is to print out some results such as number of
iteration, error value, distances and a mean distance.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 |
def golden_search(xl,xu,mode,et): it=0 e=1 while e>=et: new_x=update_interior(xl,xu) x1=new_x[0] x2=new_x[1] fx1=func_fx(m,b,x1) fx2=func_fx(m,b,x2) label=check_pos(x1,x2) #SELECTING AND UPDATING BOUNDARY-INTERIOR POINTS if mode=='max': new_boundary=find_max(xl,xu,x1,x2,label) elif mode=='min': new_boundary=find_min(xl,xu,x1,x2,label) else: print('Please define min/max mode') break #exit if mode not min or max xl=new_boundary[0] xu=new_boundary[1] xopt=new_boundary[2] distance_1=new_boundary[3] distance_2=new_boundary[4] mean_distance=(distance_1+distance_2)/2 #PLOTTING clear_output(wait=True) plt.plot([x_point,x1],[y_point,fx1],'r') plt.plot([x_point,x2],[y_point,fx2],'r') plt.plot(x,y) plt.plot(x_point,y_point,'bo') plt.axis('equal') plt.show() it+=1 print ('Iteration: ',it) r=(np.sqrt(5)-1)/2 #GOLDEN RATIO e=((1-r)*(abs((xu-xl)/xopt)))*100 #Error print('Error:',e) print ('distance 1=',distance_1) print ('distance 2=',distance_2) print('mean distance=',round(mean_distance,3)) time.sleep(1) |
We already created all functions that needed for this application. Now let's run it with the following code.
#POINT COORDINATE x_point=20 y_point=15 #START AND END LINE COORDINATES x1,y1=0,10 x2,y2=10,30 #COMPUTE m and b FOR LINE EQUATION line=line_func(x1,y1,x2,y2) m=line[0] b=line[1] x=(x1,x2) y=(y1,y2) #EXECUTE golden_search FUNCTION golden_search(0,20,'min',0.05)
If you want to try it by yourself, copy all the codes into a jupyter notebook.
If everything is fine, you should get the output as in figure 1. Anyway, I
provide the code in Jupyter notebook based on request. If you want to get it
make the
request here.