Skip to content

Commit

Permalink
Merge pull request #44 from Amir-Shamsi/numerical_analysis
Browse files Browse the repository at this point in the history
Numerical analysis
  • Loading branch information
Amir-Shamsi authored Jul 24, 2022
2 parents 56975a6 + dc8aaed commit 4886048
Show file tree
Hide file tree
Showing 8 changed files with 366 additions and 3 deletions.
7 changes: 5 additions & 2 deletions SpAlgo/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,10 @@
from SpAlgo._knapsack._knapsack import Knapsack
from SpAlgo._closest_pair._closest import ClosestPair
from SpAlgo._closest_pair._point import Point
from SpAlgo.numerical_analysis.bin_to_float import F32bit, F64bit
from SpAlgo.numerical_analysis.newton_method import Newton
from SpAlgo.numerical_analysis.secant_method import Secant

__version__ = '0.3.1'
__version__ = '0.4.0'
__all__ = ['Matrix', 'MaxSubseq', 'Array', 'Knapsack', 'ClosestPair',
'Point']
'Point', 'F32bit', 'F64bit', 'Newton', 'Secant']
111 changes: 111 additions & 0 deletions SpAlgo/numerical_analysis/bin_to_float.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
import re

class F32bit:
def __init__(self, binary_sequence: str | list) -> None:
"""
A number in 32 bit single precision IEEE 754 binary floating point standard representation requires three
building elements:
* sign (it takes 1 bit, and it's either 0 for positive or 1 for negative numbers)
* exponent (8 bits)
* mantissa (23 bits)
binary_sequence: The binary sequence in shap of string of list of 0 and 1
Example
-------
>>> bin_seq = '00000100010001000100010010111010'
>>> f = F32bit(binary_sequence=bin_seq)
>>> f.get_exponent()
>>> f.get_significand()
>>> f.get_floating_point()
"""

binary_sequence = ''.join(binary_sequence) if isinstance(binary_sequence, list) else binary_sequence

try:
if not re.match("^[01]+$", binary_sequence.strip()):
raise ValueError

self.bin_seq = [int(item) for item in binary_sequence]

if len(self.bin_seq) != 32: raise IndexError(f'Binary sequence must be 32-bits but {len(self.bin_seq)}-bits given!')

except ValueError:
raise ValueError('Please make sure the binary sequence only contains 0 and 1!')

self.__calc__()

def __calc__(self) -> None:
self._exponent = self._significand = 0

for index in range(8, 0, -1):
self._exponent += self.bin_seq[index] * (2 ** (abs(index - 8)))

for index in range(30, 8, -1):
self._significand += self.bin_seq[index] * (2 ** (-1 * (abs(index - 31))))

self._floating_point = ((-1) ** self.bin_seq[0]) * (2 ** (self._exponent - 127)) * (1 + self._significand)

def get_exponent(self) -> int | float:
return self._exponent

def get_significand(self) -> int | float:
return self._significand

def get_floating_point(self) -> int | float:
return self._floating_point


class F64bit:
def __init__(self, binary_sequence: str) -> None:
"""
A number in 64 bit single precision IEEE 754 binary floating point standard representation requires three
building elements:
* sign (it takes 1 bit, and it's either 0 for positive or 1 for negative numbers)
* exponent (11 bits)
* mantissa (52 bits)
binary_sequence: The binary sequence in shap of string of list of 0 and 1
Example
-------
>>> bin_seq = '0000010001000100010001001011101011010100010001000100010010111010'
>>> f = F64bit(binary_sequence=bin_seq)
>>> f.get_exponent()
>>> f.get_significand()
>>> f.get_floating_point()
"""
binary_sequence = ''.join(binary_sequence) if isinstance(binary_sequence, list) else binary_sequence

try:
if not re.match("^[01]+$", binary_sequence.strip()):
raise ValueError

self.bin_seq = [int(item) for item in binary_sequence]

if len(self.bin_seq) != 64: raise IndexError(f'Binary sequence must be 64-bits but {len(self.bin_seq)}-bits given!')

except ValueError:
raise ValueError('Please make sure the binary sequence only contains 0 and 1!')

self.__calc__()

def __calc__(self) -> None:
self._exponent = self._significand = 0

for index in range(11, 0, -1):
self._exponent += self.bin_seq[index] * (2 ** (abs(index - 11)))

for index in range(63, 11, -1):
self._significand += self.bin_seq[index] / (2 ** (52 - abs(index - 63)))

self._floating_point = ((-1) ** self.bin_seq[0]) * (2 ** (self._exponent - 1023)) * (1 + self._significand)

def get_exponent(self) -> int | float:
return self._exponent

def get_significand(self) -> int | float:
return self._significand

def get_floating_point(self) -> int | float:
return self._floating_point
64 changes: 64 additions & 0 deletions SpAlgo/numerical_analysis/newton_method.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
from typing import Callable, Any

class Newton:
def __init__(self, function: Callable, derivative_function: Callable) -> None:
"""Approximate solution of f(x)=0 by Newton's method.
Parameters
----------
function : `function`
Function for which we are searching for a solution f(x)=0.
derivative_function : `function`
Derivative of f(x).
Examples
--------
>>> f = lambda x: x**2 - x - 1
>>> Df = lambda x: 2*x - 1
>>> N = Newton(f, Df)
>>> N.solve(1, 1e-8, 10)
Found solution after 5 iterations.
1.618033988749989
"""

self.func = function
self.de_func = derivative_function

def solve(self, x0: int | float, epsilon: int | float, max_iter: int) -> tuple[int | float | Any, int]:
"""
Parameters
----------
x0 : `number`
Initial guess for a solution f(x)=0.
epsilon : `number`
Stopping criteria is abs(f(x)) < epsilon.
max_iter : `integer`
Maximum number of iterations of Newton's method.
Returns
-------
xn : `number`
Implement Newton's method: compute the linear approximation
of f(x) at xn and find x intercept by the formula
x = xn - f(xn)/Df(xn)
Continue until abs(f(xn)) < epsilon and return xn.
If Df(xn) == 0, return None. If the number of iterations
exceeds max_iter, then return None.
"""
xn = x0

for n in range(0, max_iter):
fxn = self.func(xn)

if abs(fxn) < epsilon:
return xn, n

de_f_xn = self.de_func(xn)

if de_f_xn == 0:
raise ZeroDivisionError('Zero derivative. No solution found.')

xn = xn - fxn / de_f_xn

raise Exception('Exceeded maximum iterations. No solution found.')
64 changes: 64 additions & 0 deletions SpAlgo/numerical_analysis/secant_method.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,64 @@
from typing import Callable


class Secant:
def __init__(self, function: Callable):
"""
Approximate solution of f(x)=0 on interval [a,b] by the secant method.
Parameters
----------
function : function
The function for which we are trying to approximate a solution f(x)=0.
Examples
--------
>>> f = lambda x: x**2 - x - 1
>>> s = Secant(f)
>>> s.solve(1, 2, 5)
1.6180257510729614
"""

self.func = function

def solve(self, a: int | float, b: int | float, N) -> int | float:
"""
Parameters
----------
a, b : numbers
The interval in which to search for a solution. The function returns
None if f(a)*f(b) >= 0 since a solution is not guaranteed.
N : (positive) integer
The number of iterations to implement.
Returns
-------
m_N : number
The x intercept of the secant line on the Nth interval
m_n = a_n - f(a_n)*(b_n - a_n)/(f(b_n) - f(a_n))
The initial interval [a_0,b_0] is given by [a,b]. If f(m_n) == 0
for some intercept m_n then the function returns this solution.
If all signs of values f(a_n), f(b_n) and f(m_n) are the same at any
iterations, the secant method fails and return None.
"""
if self.func(a) * self.func(b) >= 0:
raise Exception('Secant method fails!')

a_n, b_n = a, b

for n in range(1, N + 1):
m_n = a_n - self.func(a_n) * (b_n - a_n) / (self.func(b_n) - self.func(a_n))

f_m_n = self.func(m_n)

if self.func(a_n) * f_m_n < 0: a_n, b_n = a_n, m_n

elif self.func(b_n) * f_m_n < 0: a_n, b_n = m_n, b_n

elif f_m_n == 0: return m_n

else: raise Exception('Secant method fails!')

return a_n - self.func(a_n) * (b_n - a_n) / (self.func(b_n) - self.func(a_n))
2 changes: 1 addition & 1 deletion docs/VERSION
Original file line number Diff line number Diff line change
@@ -1 +1 @@
0.3.1
0.4.0
40 changes: 40 additions & 0 deletions src/examples/Binary to Floating Point.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,40 @@
from SpAlgo import F32bit, F64bit

# this is the knapsack algorithm.

# color code to make the output look stoning
ORANGE_COLOR = '\033[94m'

# code of make the text bold
BOLD_START = '\033[1m'

# code of making text to go back to the normal
NORM = '\033[0m'

# example of binary sequence of 32bits and 64bits
bin_seq_32bit = '00000100010001000100010010111010'
bin_seq_64bit = '0000010001000100010000001000100010001000100101110100010010111010'
"""
passing information to the class "F32bit" or "F64bit":
* binary_sequence : The function for which we are trying to approximate a solution f(x)=0.
"""
f32bit = F32bit(bin_seq_32bit)
f64bit = F64bit(bin_seq_64bit)

""" get floating points as outputs"""
float_32 = f32bit.get_floating_point()
float_64 = f64bit.get_floating_point()

"""
you can also get the values below:
* get_exponent(): returns calculated exponent
* get_significand(): returns calculated significand
"""

# we print to see the result
print('--------- F32bit & F64bit Algorithm ---------')


# printing the result
print('Calculated floating point of 32-bits binary sequence is' + BOLD_START + ORANGE_COLOR, float_32, NORM)
print('Calculated floating point of 64-bits binary sequence is' + BOLD_START + ORANGE_COLOR, float_64, NORM)
43 changes: 43 additions & 0 deletions src/examples/Newton's Method.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,43 @@
from SpAlgo import Newton

# this is the knapsack algorithm.

# color code to make the output look stoning
ORANGE_COLOR = '\033[94m'

# code of make the text bold
BOLD_START = '\033[1m'

# code of making text to go back to the normal
NORM = '\033[0m'

# example of function and its derivative
function = lambda x: x**3 - x**2 - 1
der_function = lambda x: 3*x**2 - 2*x


"""
passing information to the "Newton" class:
* function : Function for which we are searching for a solution f(x)=0.
* derivative_function : Derivative of f(x).
"""
newton = Newton(function, der_function)

""" Let's test our function newton on the polynomial to approximate the super golden ratio. """
approx, iter_count = newton.solve(x0=1, epsilon=1e-10, max_iter=10)

"""
Return Value
------------
the return value is a tuple of:
* first value: Found solution
* second value: Iteration count of found solution
"""

# we print to see the result
print('--------- Newton Algorithm ---------')


# printing the result
print('Found solution after' + BOLD_START + ORANGE_COLOR,
iter_count, NORM + 'iterations is' + BOLD_START + ORANGE_COLOR, approx, NORM)
38 changes: 38 additions & 0 deletions src/examples/Secant Method.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,38 @@
from SpAlgo import Secant

# this is the knapsack algorithm.

# color code to make the output look stoning
ORANGE_COLOR = '\033[94m'

# code of make the text bold
BOLD_START = '\033[1m'

# code of making text to go back to the normal
NORM = '\033[0m'

# example of function and its derivative
polynomial = lambda x: x**3 - x**2 - 1

"""
passing information to the "Secant" class:
* function : The function for which we are trying to approximate a solution f(x)=0.
"""
secant = Secant(polynomial)

""" Let's test our function with input values for which we know the correct output. Let's find an approximation of the
super golden ratio: the only real root of the polynomial """
approx = secant.solve(a=1, b=2, N=20)

"""
Return Value
------------
the return value is Found solution
"""

# we print to see the result
print('--------- Secant Algorithm ---------')


# printing the result
print('Found solution is' + BOLD_START + ORANGE_COLOR, approx, NORM)

0 comments on commit 4886048

Please sign in to comment.