num-an

num-an

Numerical Analysis Methods Implementations

.gitignore

__pycache__/
test.py

README.md

To be done:

• Documentation
• Bugfixes

readme.txt

verbose:
0: print nothing
1: print enter, exit messages, O(1) computing steps
2: print O(n) computing steps
3: print All computing steps
4: print all info and formulas

Interactive:
False: disabled
True: stop after every step

Messages:
~~~Operation being performed~~~
###Variable name###
***Expected action***

solve_library.py

"""This module ilustrates some methods of numerical analysis"""
__author__ = "Oleg Zubchenko"

from math import *
import numpy
import pdb
from copy import deepcopy
from functools import reduce
from itertools import chain
import inspect

numpy.set_printoptions(precision = 3, suppress = True)

def pause():
	"""Waits for enter key press"""
	input("***Hit [Enter]***")

def fun_header():
	"""Prints function start message"""
	print("~~~Begin %s~~~" % inspect.stack()[1][3])

def fun_footer():
	"""Prints function end message"""
	print("~~~End %s~~~" % inspect.stack()[1][3])

def print_group(*Matrices):
	"""Prints matrices and vectors side by side"""
	if len(Matrices) == 0: return
	for M in Matrices:
		assert(M.shape[0] == Matrices[0].shape[0])
	order = Matrices[0].shape[0]
	M_lines = []
	for M in Matrices:
		M_lines.append(str(M)[1:-1].split('\n '))
	lines = [""] * order
	for i in range(order):
		lines[i] = "(" + ("|".join([
			M_lines[j][i][1:-1] for j in range(len(Matrices))
		])) + ")"
	print('\n'.join(lines))

def swap_rows(M, i, j):
	"""Swaps matrix rows if possible"""
	if not i in range(M.shape[0]):
		pdb.set_trace()
		raise IndexError('index %d out of bounds for axis 0 with size %d' % (i, n))
	if not j in range(M.shape[0]):
		pdb.set_trace()
		raise IndexError('index %d out of bounds for axis 0 with size %d' % (j, n))
	M[[i, j], :] = M[[j, i], :]

def swap_cols(M, i, j):
	"""Swaps matrix cols if possible"""
	if not i in range(M.shape[1]):
		pdb.set_trace()
		raise IndexError('index %d out of bounds for axis 0 with size %d' % (i, n))
	if not j in range(M.shape[1]):
		pdb.set_trace()
		raise IndexError('index %d out of bounds for axis 0 with size %d' % (j, n))
	M[:, [i, j]] = M[:, [j, i]]

def to_triangle_top(A, *, do_swap_rows = False, do_swap_cols = False,
		verbose = 0, interactive = False):
	"""Transforms given matrix to upper triangle form
	Algorithm:
		PA = U
		PAQ = U
	Args:
		A (array like) - matrix to be transformed
		no_swap (bool) - do not swap cols
	Returns:
		(P, Q)
	"""
	A = numpy.array(A).astype(float)
	assert (A.shape[0] == A.shape[1])
	if verbose >= 1:
		fun_header()
		print("###A###"), print_group(A)
		if interactive: pause()
	n = A.shape[0]
	P = numpy.eye(n).astype(float)
	Q = numpy.eye(n).astype(float)
	for i in range(n):
		if do_swap_rows and do_swap_cols:
			swapee = numpy.unravel_index(abs(A)[i:, i:].argmax(), (n-i, n-i))
			swapee = tuple([i + idx for idx in swapee])
			if abs(A).item(swapee) < 1e-9:
				if verbose >= 1:
					print("~~~Singular matrix~~~")
					print("###P|A|Q###"),print_group(P, A, Q)
					if interactive: pause()
				raise ValueError("singular matrix")
			swap_rows(P, i, swapee[0])
			swap_rows(A, i, swapee[0])
			if verbose >= 2:
				print("~~~swap rows %d %d~~~" % (i, swapee[0]))
				print("###P|A|Q###"),print_group(P, A, Q)
				if interactive: pause()
			swap_cols(Q, i, swapee[1])
			swap_cols(A, i, swapee[1])
			if verbose >= 2:
				print("~~~swap cols %d %d~~~" % (i, swapee[1]))
				print("###P|A|Q###"),print_group(P, A, Q)
				if interactive: pause()
		elif do_swap_rows:
			swapee = numpy.unravel_index(abs(A)[i:, i].argmax(), (n-i, 1))
			swapee = tuple([i + idx for idx in swapee])
			swap_rows(P, i, swapee[0])
			swap_rows(A, i, swapee[0])
			if verbose >= 2:
				print("~~~swap rows %d %d~~~" % (i, swapee[0]))
				print("###P|A|Q###"),print_group(P, A, Q)
				if interactive: pause()
		elif do_swap_cols:
			swapee = numpy.unravel_index(abs(A)[i, i:].argmax(), (1, n-i))
			swapee = tuple([i + idx for idx in swapee])
			swap_cols(Q, i, swapee[1])
			swap_cols(A, i, swapee[1])
			if verbose >= 2:
				print("~~~swap cols %d %d~~~" % (i, swapee[1]))
				print("###P|A|Q###"),print_group(P, A, Q)
				if interactive: pause()
		for j in range(i + 1, n):
			#pdb.set_trace()
			P[j] -= P[i] * A[j][i] / A[i][i]
			A[j] -= A[i] * A[j][i] / A[i][i]
			if verbose >= 3:
				print("~~~Substracting row %d from row %d~~~" % (i, j))
				print("###P|A|Q###"), print_group(P, A, Q)
				if interactive: pause()
		if verbose >= 2:
			print("~~~Substracting row %d from the below~~~" % i)
			print("###P|A|Q###"), print_group(P, A, Q)
			if interactive: pause()
	if verbose >= 1:
		print("~~~Finished transform to top triangle~~~")
		print("###P|A|Q###"),print_group(P, A, Q)
		if interactive: pause()
		fun_footer()
	return P, Q

#(A_triangle) * (Q^-1 * x) = (b) => (x)
def solve_triangle_top(U, Q, b, *, verbose = 0, interactive = False):
	"""Solves upper triangle system
	Algorithm:
		(U) (Q^-1 x) = f
		Q^-1 x = y
		U y = f -> find y
		x = Q y
	"""
	U = numpy.array(U).astype(float)
	if Q == None:
		Q = numpy.identity(U.shape[0])
	Q = numpy.array(Q).astype(float)
	b = numpy.array(b).astype(float)
	assert (Q.shape[0] == b.shape[0])
	assert (b.shape[1] == 1)
	if verbose >= 1:
		fun_header()
		print("###U|Q|b###")
		print_group(U, Q, b)
		if interactive: pause()
	n = U.shape[0]
	if verbose >= 1:
		print("~~~Composing system~~~")
		print("###U|b###")
		print_group(U, b)
		if interactive: pause()
	for i in reversed(range(n)):
		b[i] /= U[i, i]
		U[i] /= U[i, i]
		if verbose >= 2:
			print("~~~Normalizing row %d~~~" % i)
			print_group(U, b)
			if interactive: pause()
		for j in range(i):
			b[j] -= b[i] * U[j, i]
			U[j] -= U[i] * U[j, i]
			if verbose >= 3:
				print("~~~Substracting row %d from %d~~~" % (i, j))
				print_group(U, b)
				if interactive: pause()
		if verbose >= 2:
			print("~~~Substracting row %d from the above~~~" % i)
			print_group(U, b)
			if interactive: pause()
	if verbose >= 1:
		print("~~~System solved~~~")
		print("###U|Q^-1|b###")
		print_group(U, Q.T, b)
		if interactive: pause()
	x = numpy.dot(Q, b)
	if verbose >= 1:
		print("~~~Reversing column permutation~~~")
		print("###x###")
		print_group(x)
		if interactive: pause()
		fun_footer()
	return x

def solve_triangle_bottom(L, Q, b, *, verbose = 0, interactive = False):
	"""Solves bottom triangle system
	Algorithm:
		(L) (Q^-1 x) = f
		Q^-1 x = y
		L y = f -> find y
		x = Q y
	"""
	L = numpy.array(L).astype(float)
	if Q == None:
		Q = numpy.identity(L.shape[0])
	Q = numpy.array(Q).astype(float)
	b = numpy.array(b).astype(float)
	assert (L.shape == Q.shape)
	assert (Q.shape[0] == b.shape[0])
	#assert (b.shape[1] == 1)
	if verbose >= 1:
		fun_header()
		print("###L|Q|b###")
		print_group(L, Q, b)
		if interactive: pause()
	n = L.shape[0]
	if verbose >= 1:
		print("~~~Composing system~~~")
		print("###L|b###")
		print_group(L, b)
		if interactive: pause()
	for i in range(n):
		b[i] /= L[i, i]
		L[i] /= L[i, i]
		if verbose >= 2:
			print("~~~Normalizing row %d~~~" % i)
			print_group(L, b),
			if interactive: pause()
		for j in range(i+1, n):
			b[j] -= b[i] * L[j, i]
			L[j] -= L[i] * L[j, i]
			#if verbose >= 3:
			#	print_group(L, b)
			#	print("~~~Substracting row %d from %d~~~" % (i, j))
			#	if interactive: pause()
		if verbose >= 2:
			print("~~~Substracting row %d from the below~~~" % i)
			print_group(L, b)
			if interactive: pause()
	if verbose >= 1:
		print("~~~System solved~~~")
		print("###Q^-1|b###")
		print_group(Q.T, b)
	x = numpy.dot(Q, b)
	if verbose >= 1:
		print("~~~Reverting column permutation~~~")
		print("###x###")
		print_group(x)
		if interactive: pause()
		fun_footer()
	return x

def solve_gauss(A, b, *, do_swap_rows = False, do_swap_cols = False,
		verbose = 0, interactive = False):
	"""Solves system using gaussian ellimination
	Algorithm:
		Forward elimination:
			P A Q = U
		Backward elimination:
			U (Q^-1 x) = (P b)
			(Q^-1 x) = y
		Reversing permutation:
		x = Q y
	"""
	A = numpy.array(A).astype(float)
	b = numpy.array(b).astype(float)
	assert (A.shape[0] == A.shape[1])
	assert (A.shape[0] == b.shape[0])
	assert (b.shape[1] == 1)
	if verbose >= 1:
		fun_header()
		print("###A|b###")
		print_group(A, b)
	P, Q = to_triangle_top(A, do_swap_rows=do_swap_rows,
			do_swap_cols=do_swap_cols, verbose = verbose, interactive = interactive)
	if verbose >= 1:
		print("~~~Done forward elimination~~~")
		U = numpy.dot(numpy.dot(P, A), Q)
		b = numpy.dot(P, b)
		print("###U|b###")
		print_group(U, b)
	x = solve_triangle_top(U, Q, b, verbose = verbose, interactive = interactive)
	if verbose >= 1:
		print("~~~Done backward elimination~~~")
		print("###x###")
		print_group(x)
		fun_footer()
	return x

def solve_gauss_noswap(A, b, *, verbose = 0, interactive = False):
	return solve_gauss(A, b, do_swap_rows = False, do_swap_cols = False,
			verbose = verbose, interactive = interactive)

#Archived
#def lu_decompose_noswap(A, *, verbose = 0, interactive = False):
#	assert (A.shape[0] == A.shape[1])
#	order = A.shape[0]
#	L = numpy.identity(order)
#	if verbose >= 1: print_group(L, A), print("~~~Started  LU decomposition~~~")
#	for i in range(order):
#		if verbose >= 2: print_group(L, A), print("~~~Building L col %d~~~" % i)
#		for j in range(i + 1, order):
#			if verbose >= 3:
#				print_group(L, A), print("~~~Building L[%d, %d]~~~" % (j, i))
#			L[j, i] = A[j, i] / A[i, i]
#			A[j] -= A[i] / A[i, i] * A[j, i]
#	if verbose >= 1: print_group(L, A), print("~~~Finished LU decomposition~~~")
#	return L

def lu_decompose(A, *, verbose = 0, interactive = False):
	"""Does LU decomposition for given matrix:
	Algorithm: Actually, gaussian elimination with memorising each step
		L_n-1 P_n-1 ... L_0 P_0 A Q_0 ... Q_n-1 = U
		L_n-1 P_n-1 ... L_0 P_0 
                    = L'_n-1 ... L'_0 P'_n-1 ... P'_0 = L_composed P
		L_composed P A Q = U // L_composed to rhs
		P A Q = L U
	"""
	A = numpy.array(A).astype(float)
	assert (A.shape[0] == A.shape[1])
	n = A.shape[0]
	Ls = []
	Ps = []
	Qs = []
	U = A
	if verbose >= 1:
		fun_header()
		print("###A###")
		print_group(A)
		if interactive: pause()
	for i in range(n-1):
		if verbose >= 3:
			print("~~~Rounding error: %f~~~" % numpy.linalg.norm(
					reduce(numpy.dot,
					chain(chain(*zip(reversed(Ls), reversed(Ps))), [A], Qs),
					numpy.identity(n)) - U))
			print("###Ls###"), print_group(*reversed(Ls))
			print("###Ps###"), print_group(*reversed(Ps))
			print("###A###"), print_group(A)
			print("###Qs###"), print_group(*Qs)
			if interactive: pause()
		if verbose >= 2:
			print("###U###"), print_group(U)
			print("~~~Elliminating col %d~~~" % i)
			if interactive: pause()
		swapee = numpy.unravel_index(abs(U)[i:, i:].argmax(), (n-i, n-i))
		swapee = tuple([i + idx for idx in swapee])
		if (abs(U.item(swapee)) < 1e-9):
			if verbose >= 1:
				print("~~~Singular matrix~~~")
				pdb.set_trace()
			raise ValueError("Singular matrix")
		Ls.append(numpy.identity(n))
		Ps.append(numpy.identity(n))
		Qs.append(numpy.identity(n))
		swap_rows(Ps[i], i, swapee[0])
		swap_cols(Qs[i], i, swapee[1])
		if verbose >= 2:
			print("~~~Swapping cols %d with %d~~~" % (i, swapee[1]))
			print("###U, Qs[%d]###" % i)
			print_group(U, Qs[i])
			if interactive: pause()
		U = numpy.dot(U, Qs[i])
		if verbose >= 2:
			print("~~~Swapping rows %d with %d~~~" % (i, swapee[0]))
			print("###Ps[%d], U###" % i)
			print_group(Ps[i], U)
			if interactive: pause()
		U = numpy.dot(Ps[i], U)
		for j in range(i+1, n):
			Ls[i][j, i] = -U[j, i] / U[i, i]
		if verbose >= 2:
			print("~~~Zeroing col %d of U~~~" % i)
			print("###Ls[%d], U###" % i)
			print_group(Ls[i], U)
			if interactive: pause()
		U = numpy.dot(Ls[i], U)
	if verbose >= 2:
		print("~~~Done transforming~~~")
		print("###Ls###"), print_group(*reversed(Ls))
		print("###Ps###"), print_group(*reversed(Ps))
		print("###A###"), print_group(A)
		print("###Qs###"), print_group(*Qs)
		if interactive: pause()
	if verbose >= 1:
		print("~~~Rearranging L and P matrices~~~")
		print("###U###"), print_group(U)
		if interactive: pause()
	L_bars = []
	for i in range(n-1):
		Tangle = reduce(numpy.dot, Ps[:i:-1], numpy.identity(n))
		L_bars.append(numpy.dot(numpy.dot(Tangle, Ls[i]), Tangle.T))
		if verbose >= 2:
			print("###Left permutation matrix for L_bar[%d]###" % i)
			print_group(*Ps[:i:-1])
			print("###Left/right permutation composed [%d]###" % i)
			print_group(Tangle, Tangle.T)
			print("###L_bar[%d]###" %i)
			print_group(L_bars[i])
			if interactive: pause()
	P = reduce(numpy.dot, reversed(Ps), numpy.identity(n))
	Q = reduce(numpy.dot, Qs, numpy.identity(n))
	L_rhs = numpy.linalg.inv(
			reduce(numpy.dot, reversed(L_bars), numpy.identity(n))
	)
	if verbose >= 1:
		print("~~~Composing Resulting matrices~~~")
		print("###P|A|Q=L|U###"), print_group(P, A, Q, L_rhs, U)
		print("Finished LU full swap decomposition~~~")
		print("~~~Rounding error: %f~~~" % numpy.linalg.norm(
				numpy.dot(numpy.dot(P, A), Q) - numpy.dot(L_rhs, U)
		))
		if interactive: pause()
	return (P, Q, L_rhs, U)

def solve_lu(A, b, *, verbose = 0, interactive = False):
	"""Solves system using LU decomposition.
	If you only need to decompose matrix, see lu_decompose
	Algorithm:
		A * x = b
		(P * A * Q) * (Q^-1 * x) = P * b
		L * U * (Q^-1 * x) = P * b
		L * (U * (Q^-1 * x)) = (P * b)
		U * (Q^-1 * x) = c = <solve bottom triangle(L, None, P * b)>
		U * Q^-1 * x = c
		x = <solve top triangle(U, Q, c)>
	"""
	A = numpy.array(A).astype(float)
	b = numpy.array(b).astype(float)
	assert (A.shape[0] == A.shape[1])
	assert (A.shape[0] == b.shape[0])
	assert (b.shape[1] == 1)
	if verbose >= 1:
		fun_header()
		print("###A|b###")
		print_group(A, b)
		if interactive: pause()
	n = A.shape[0]
	P, Q, L, U = lu_decompose(A, verbose = verbose, interactive = interactive)
	c = solve_triangle_bottom(L, None, numpy.dot(P, b), verbose = verbose,
			interactive = interactive)
	if verbose >= 1:
		print("~~~Done solving L system~~~")
		print("~~~Solving U system~~~")
		if interactive: pause()
	x = solve_triangle_top(U, Q, c, verbose = verbose, interactive = interactive)
	if verbose >= 1:
		print("~~~Finished solving using COMPLETE LU decomposition~~~")
		print("###x###")
		print(x)
		if interactive: pause()
		fun_footer()
	return x

def sqrt_decompose(A, *, verbose = 0, interactive = False):
	"""Does sqrt decomposition for a given matrix
	Find L: A = L L^T
	Requirements:
		A = A^T
		A >> 0 - positive definite matrix
	"""
	A = numpy.array(A).astype(float)
	assert (A.shape[0] == A.shape[1])
	n = A.shape[0]
	if verbose >= 1:
		fun_header()
		print("###A###")
		print_group(A)
		if interactive: pause()
	L = numpy.identity(n)
	for i in range(n):
		L[i,i] = sqrt(A[i, i] - sum(L[i, :i] * L[i, :i]))
		if verbose >= 2:
			print("###L[%d, %d] = sqrt(%.3f - (%s)) = %.3f" % (i, i, A[i, i],
					" + ".join(["%.3f" % k for k in L[i, :i] * L[i, :i]]),
					L[i, i]
			))
			if interactive: pause()
		for j in range(i+1, n):
			L[j, i] = (A[j, i] - sum(L[i, :i] * L[j, :i])) / L[i, i]
			if verbose >= 2:
				print("###L[%d, %d] = (%.3f - (%s)) / %.3f = %.3f" % (j, i, A[j, i],
						" + ".join(["%.3f" % k for k in L[i, :i] * L[j, :i]]),
						L[i, i], L[j, i]
				))
				if interactive: pause()
	if verbose >= 1:
		print("~~~Finished sqrt decomposition~~~")
		print("###L*L^T=A###")
		print_group(L, L.T, numpy.dot(L, L.T))
		if interactive: pause()
		fun_footer()
	return L

def solve_sqrt(A, b, *, verbose = 0, interactive = False):
	"""Solves sytem using sqrt decomposition
	Algorithm:
		A x = b
		L L^T x = b
		L^T x = y
		solve L y = b, get y
		solve L^T x = y, get x
	"""
	A = numpy.array(A).astype(float)
	b = numpy.array(b).astype(float)
	assert (A.shape[0] == A.shape[1])
	assert (A.shape[0] == b.shape[0])
	assert (b.shape[1] == 1)
	if verbose >= 1:
		fun_header()
		print("###A|b###")
		print_group(A, b)
		if interactive: pause()
	n = A.shape[0]
	L = sqrt_decompose(A, verbose = verbose, interactive = interactive)
	y = solve_triangle_bottom(L, None, b, verbose=verbose,
			interactive=interactive)
	if verbose >= 1:
		print("~~~Finished solve_sqrt step 1(bottom)~~~")
		print("###L * y = b###")
		print_group(L, y, b)
		if interactive: pause()
	x = solve_triangle_top(L.T, None, y, verbose=verbose, interactive=interactive)
	if verbose >= 1:
		print("~~~Finished solve_sqrt step 2(top)~~~")
		print_group(A, x, b)
		if interactive: pause()
		fun_footer()
	return x

def solve_tridiagonal(A, b, *, verbose = 0, interactive = False):
	"""Solves tridiagonal system using tridiagonal method
	Algorithm:
		c[0, n-1] = diag(A)
		b[0, n-2] = -diag(A, +1)
		a[1, n-1] = -diag(A, -1)
		f[0, n-1] = flatten(b)
		alpha[1] = b[0] / c[0]
		beta[1] = f[0] / c[0]
		alpha[i+1] = b[i] / (c[i] - a[i] alpha[i])
		beta[i+1] = (f[i] + a[i] beta[i]) / (c[i] - a[i] alpha[i])
		to alpha[n-1], beta[n-1]
		x[n-1] = (f[n-1] + a[n-1] beta[n-1]) / (c[n-1] - alpha[n-1] a[n-1])
		x[i] = alpha[i+1] x[i+1] + beta[i+1]
		to x[0]
	"""
	A = numpy.array(A).astype(float)
	b = numpy.array(b).astype(float)
	assert (A.shape[0] == A.shape[1])
	assert (A.shape[0] == b.shape[0])
	assert (b.shape[1] == 1)
	if verbose >= 1:
		fun_header()
		print("###A|b###")
		print_group(A, b)
		if interactive: pause()
	f = b.flatten()
	n = A.shape[0]
	c = numpy.diag(A, 0)
	b = -numpy.diag(A, 1)
	a = -numpy.concatenate(([0.], numpy.diag(A, -1)) )
	alpha = numpy.zeros(n)
	beta  = numpy.zeros(n)
	alpha[1] = b[0] / c[0]
	beta [1] = f[0] / c[0]
	if verbose >= 2:
		print("alpha[%d] = %f / %f = %f" % (1, b[0], c[0], alpha[1]))
		print("beta [%d] = %f / %f = %f" % (1, b[0], c[0], alpha[1]))
	for i in range(1, n-1):
		alpha[i+1] =                    b[i] / (c[i] - a[i] * alpha[i])
		beta [i+1] = (f[i] + a[i] * beta[i]) / (c[i] - a[i] * alpha[i])
		if verbose >= 2:
			print("alpha[%d] = %f / (%f - %f * %f) = %f" % 
					(i+1, b[i], c[i], a[i], alpha[i], alpha[i+1])
			)
			print("beta [%d] = (%f + %f * %f) / (%f - %f * %f) = %f" % 
					(i+1, f[i], a[i], b[i], c[i], a[i], alpha[i], beta[i+1])
			)
			if interactive: pause()
	x = numpy.zeros(n)
	x[n-1] = (f[n-1] + a[n-1] * beta[n-1]) / (c[n-1] - alpha[n-1] * a[n-1])
	if verbose >= 2:
		print("x[%d] = (%f + %f * %f) / (%f - %f * %f) = %f" %
				(n-1, f[n-1] ,a[n-1], beta[n-1], c[n-1], alpha[n-1], a[n-1], x[n-1])
		)
		if interactive: pause()
	for i in reversed(range(n-1)):
		x[i] = alpha[i+1] * x[i+1] + beta[i+1]
		if verbose >= 2:
			print("x[%d] = %f * %f + %f = %f" %
					(i, alpha[i+1], x[i+1], beta[i+1], x[i])
			)
			if interactive: pause()
	x = x.reshape(n, 1)
	if verbose >= 2:
		print("~~~System solved~~~")
		print("###x###")
		print_group(x)
		fun_footer()
	return x

def solve_jacobi(A, b, *, eps=1e-3, weight=1., verbose=0, interactive=False):
	"""Solves system with given precision by Jacobi method
	Algorithm:
	A * x = b
	D * A_unite * x = b
	A_unite * x = D^-1 * b
	(A_unite - E + E) * x = D^-1 * b
	E x = (E - A_unite) * x + D^-1 * b
		E - A_unite = B
		D^-1 b = g
	x = B * x + g
	Requirements:
		B - a contraction mapping
	"""
	A = numpy.array(A).astype(float)
	b = numpy.array(b).astype(float)
	assert (A.shape[0] == A.shape[1])
	assert (A.shape[0] == b.shape[0])
	assert (b.shape[1] == 1)
	n = A.shape[0]
	if verbose >= 1: print("~~~Started solve using jacobi method~~~")
	D = numpy.diag(A).reshape(n, 1)
	g = b / D#::= g
	B = numpy.identity(n) - A / D
	q = numpy.linalg.norm(B)
	target_eps = abs(eps) * (1 - q)
	weight = 0. if weight <= 0. else (1. if weight >= 1. else float(weight))
	if verbose >= 1: print("###q = %f|eps = %f|weight = %f###" %
			(q, target_eps, weight)
	)
	if (q >= 1):
		if verbose >= 1: print("!!!NON CONVERGENT SEQUENCE!!!")
		return None
	if verbose >= 1:
		print("###B|g###"), print_group(B, g)
		print("~~~Found iterative metod for A~~~")
		if interactive: pause()
	x = deepcopy(g)
	k = 0
	while True:
		x_next = weight * (numpy.dot(B, x) + g) + (1 - weight) * x
		eps = numpy.amax(numpy.absolute(x_next - x))
		if verbose >= 2:
			print_group(B, x, g, x_next), print("###B, x_%d, g, x_%d###" % (k, k + 1))
			print("###eps = %f###" % eps)
			print("~~~Applying iteration %d~~~" % k)
			if interactive: pause()
		k += 1
		x = x_next
		if eps < target_eps:
			if verbose >= 1:
				print("###x###"), print(x)
				print("~~~Finished solve jacobi method~~~")
				return x
	print("???Unreachable code solve_jacobi???")

def solve_seidel(A, b, *, eps = 1e-3, verbose = 0, interactive = False):
#	[A] * x = [b]
#	D * A_unite * x = b
#	A_unite * x = D^-1 * b
#	(E - E + A) * x = D^-1 * b
#	x = (A - E) * x + D^-1 * b
#	x = B * x + g
#	for every coordinate separately
	A = deepcopy(A)
	b = deepcopy(b)
	assert (A.shape[0] == A.shape[1])
	assert (A.shape[0] == b.shape[0])
	assert (b.shape[1] == 1)
	order = A.shape[0]
	if verbose >= 1: print("~~~Started solve using zeidel method~~~")
	D = numpy.diag(A).reshape(order, 1)
	g = b / D
	B = numpy.identity(order) - A / D
	q = numpy.linalg.norm(B)
	target_eps = abs(eps) * (1 - q)
	if verbose >= 1: print("##q = %f|eps = %f###" % (q, target_eps))
	if (q >= 1):
		if verbose >= 1: print("!!!NON CONVERGENT SEQUENCE!!!")
		return None
	if verbose >= 1:
		print("##B|g###"), print_group(B, g),
		print("~~~Found iterative metod for A~~~")
		if interactive: pause()
	x = deepcopy(g)
	i = 0
	while True:
		i += 1
		if verbose >= 2: print("~~~Building approx %d~~~" % i)
		x_prev = deepcopy(x)
		for j in range(order):
			if verbose >= 2:
				print("x_%d[%d] = %s + %f = %f" % (i, j,
						' '.join(["%+.3f" % i for i in B[j] * x.reshape(4)]),
						g[j], numpy.dot(B[j], x) + g[j]))
			x[j] = numpy.dot(B[j], x) + g[j]
		eps = numpy.amax(numpy.absolute(x - x_prev))
		if verbose >= 2:
			print_group(B, x_prev, g, x)
			print("###x###"), print(x)
			print("###eps = %f###" % eps)
			print("~~~Done iteration %d~~~" % i)
			if interactive: pause()
		if eps < target_eps:
			if verbose >= 1:
				print("###x###"), print(x)
				print("~~~Finished solve seidel method~~~")
				return x
	print("???Unreachable code solve_zeidel???")

def solve_gauss_seidel(A, b, *, eps = 1e-3, verbose = 0, interactive = False):
	#[A] * x = [b]
	#(L + U) * x = b
	#L * x = b - U * x
	#x = L^-1(b - U * x)
	#x = (-L^-1 * U) * x + (L^-1 * b)
	A = deepcopy(A)
	b = deepcopy(b)
	assert (A.shape[0] == A.shape[1])
	assert (A.shape[0] == b.shape[0])
	assert (b.shape[1] == 1)
	order = A.shape[0]
	if verbose >= 0:
		print_group(A, b), print("~~~Started solving system with gauss-seidel~~~")
	L = numpy.tril(A)
	U = numpy.triu(A, 1)
	L_inv = solve_triangle_bottom(L, None, numpy.identity(order), verbose=verbose,
			interactive = interactive)
	B = -numpy.dot(L_inv, U)
	g = numpy.dot(L_inv, b)
	q = numpy.linalg.norm(B)
	target_eps = abs(eps) * (1 - q)
	if verbose >= 1: print("##q = %f|eps = %f###" % (q, target_eps))
	if (q >= 1):
		if verbose >= 1: print("!!!NON CONVERGENT SEQUENCE!!!")
		return None
	if verbose >= 1:
		print("##B|g###"), print_group(B, g)
		print("~~~Found iterative metod for A~~~")
		if interactive: pause()
	x = deepcopy(g)
	i = 0
	while True:
		i += 1
		if verbose >= 2: print("~~~Building approx %d~~~" % i)
		x_prev = deepcopy(x)
		for j in range(order):
			if verbose >= 2:
				print("x_%d[%d] = %s + %f = %f" % (i, j,
						' '.join(["%+.3f" % i for i in B[j] * x.reshape(4)]),
						g[j], numpy.dot(B[j], x) + g[j]))
			x[j] = numpy.dot(B[j], x) + g[j]
		eps = numpy.amax(numpy.absolute(x - x_prev))
		if verbose >= 2:
			print_group(B, x_prev, g, x)
			print("###x###"), print(x)
			print("###eps = %f###" % eps)
			print("~~~Done iteration %d~~~" % i)
			if interactive: pause()
		if eps < target_eps:
			if verbose >= 1:
				print("###x###"), print(x)
				print("~~~Finished solve gauss seidel method~~~")
				return x
	print("???Unreachable code solve_gauss_seidel???")

def solve_grad_decent(A, b, *, eps=1e-3, verbose=0, interactive=False,
		simple_guess=False):
	A = deepcopy(A)
	b = deepcopy(b)
	assert (A.shape[0] == A.shape[1])
	assert (A.shape[0] == b.shape[0])
	assert (b.shape[1] == 1)
	order = A.shape[0]
	if verbose >= 0: 
		print_group(A, b), print("~~~Started solving system with grad-decent~~~")
	if simple_guess:
		x = numpy.zeros((order, 1)).astype(float)
		x[0, 0] = 1.
	else:
		x = numpy.random.random_sample((order, 1))
		x /= numpy.linalg.norm(x)
	target_eps = eps
	k = 0
	while True:
		r = numpy.dot(A, x) - b
		x_next = x - numpy.dot(r.T, r) / numpy.dot(numpy.dot(A, r).T, r) * r
		eps = numpy.linalg.norm(x - x_next)
		if verbose >= 2:
			print_group(A, x, b, r, x_next), print("###A, x, b, r, x_next###")
			print("###eps = %f###" % eps)
			print("~~~Done iteration %d~~~" % k)
			if interactive: pause()
		k += 1
		x = x_next
		if eps < target_eps:
			print(x), print("~~~Finished grad decent~~~")
			return x
	print("???Unreachable code solve_gradient_decent???")

def solve_coordinate_decent(A, b, *, eps=1e-3, verbose=0, interactive=False,
		simple_guess=False):
	A = deepcopy(A)
	b = deepcopy(b)
	assert (A.shape[0] == A.shape[1])
	assert (A.shape[0] == b.shape[0])
	assert (b.shape[1] == 1)
	order = A.shape[0]
	if verbose >= 0:
		print_group(A, b),
		print("~~~Started solving system with coordinate-decent~~~")
	if simple_guess:
		x = numpy.zeros((order, 1)).astype(float)
		x[0, 0] = 1.
	else:
		x = numpy.random.random_sample((order, 1))
		x /= numpy.linalg.norm(x)
	target_eps = eps
	i = 0
	while True:
		i += 1
		x_prev = deepcopy(x)
		for j in range(order):
			r = numpy.dot(A, x) - b
			x[j] = x[j] - numpy.dot(r.T, r) / numpy.dot(numpy.dot(A, r).T, r) * r[j]
			if verbose >=2:
				print("x_%d[%d] = %f" % (i, j, x[j]))
				if verbose >= 3:
					if interactive: pause()
		eps = numpy.linalg.norm(x - x_prev)
		if verbose >= 2:
			print("###eps = %f###" % eps)
			print("~~~Done iteration %d~~~" % i)
			if interactive: pause()
		if eps < target_eps:
			if verbose >= 2:
				print("###x###"), print(x), print("~~~Finished coordinate decent~~~")
			return x
	print("???Unreachable code solve_gradient_decent???")

def rotation_matrix(order, i, j, theta):
	t = theta
	result = numpy.identity(order)
	result[[[i], [j]], [i, j]] = numpy.array(
			[cos(t), -sin(t), sin(t), cos(t)]
	).reshape(2, 2)
	return result

def eig_rotation(A, *, eps = 1e-3, verbose = 0, interactive = False):
	A = deepcopy(A)
	target_eps = abs(eps)
	assert(A.shape[0] == A.shape[1])
	order = A.shape[0]
	eps = 2 * target_eps
	k = 0
	H = numpy.identity(order)
	if verbose >= 1:
		print(A), print("~~~Started computing eigens using rotations~~~")
	while True:#current_eps > eps
		i, j = numpy.unravel_index(numpy.triu(abs(A), 1).argmax(), A.shape)
		eps = abs(A.item(i, j))
		if eps <= target_eps:
			if verbose >= 1:
				print(A),
				print("~~~Finished computing eigens using rotations (k=%d)~~~" % k)
			return (numpy.diag(A), H)
		phi = 1 / 2 * atan(2 * A.item(i, j) / (A[i, i] - A[j, j]))
		H_step = rotation_matrix(order, i, j, phi)
		H = numpy.dot(H, H_step)
		if verbose >= 2:
			print("###Iteration %d(eps = %f)###" % (k, eps))
			print("###H_%d^T|A_%d|H_%d###" % (k, k, k))
			print_group(H_step.T, A, H_step)
			print("~~~Applying R_%d_%d(%f)~~~" % (i, j, phi))
			if interactive: pause()
		A = numpy.dot(numpy.dot(H_step.T, A), H_step)
		k += 1
	return "unreachable code eig_rotation"

# A => (x, lambda)
def eig_max_power(A, *, eps=1e-3, verbose=0, interactive=False,
		simple_guess=False):
	A = deepcopy(A)
	target_eps = abs(eps)
	assert(A.shape[0] == A.shape[1])
	order = A.shape[0]
	if verbose >= 1:
		print(A), print("~~~Started computing dominant eigen using power method~~~")
	if simple_guess:
		x = numpy.zeros((order, 1)).astype(float)
		x[0, 0] = 1.
	else:
		x = numpy.random.random_sample((order, 1))
		x /= numpy.linalg.norm(x)
	eps = float('inf')
	lambda_ = float('inf')
	deviation = float('inf')
	deviation_eps = float('inf')
	k = 0
	while True:
		x_next = numpy.dot(A, x)
		lambda_next = numpy.linalg.norm(x_next)
		deviation_next = numpy.linalg.norm(x_next / lambda_next - x)
		eps = abs(lambda_next - lambda_)
		deviation_eps = abs(deviation_next - deviation)
		if verbose >= 2:
			print("###A|x_%d|x_%d|normalized(x_%d)|deviation###" % (k, k+1, k+1))
			print_group(A, x, x_next, x_next / lambda_next, x_next / lambda_next - x)
			print("###k = %d|lambda = %f|eps = %f|dev = %f|dev_eps = %f###" %
					(k, lambda_next, eps, deviation_next, deviation_eps)
			)
			print("~~~Applying iteration %d~~~" % k)
			if interactive: pause()
		k += 1
		x = x_next / lambda_next
		lambda_ = lambda_next
		deviation = deviation_next
		if eps <= target_eps:
			if deviation > eps:
				if verbose >= 1:
					print("!!!WARNING: EIGENVECTOR IS UNTRUSTED!!!")
				return (numpy.array([lambda_next]), numpy.zeros((order, 1)))
			if verbose >= 1:
				print("###x###"), print(x_next)
				print("~~~Finished computing dominant eig using power method~~~")
			return (numpy.array([lambda_next]), x_next)
	return "unreachable code eig_power"

def polynom_krylov(A, *, solver=solve_gauss_noswap, verbose=0,
		interactive=False):
	#A
	#y_0 = [1, 0, ..., 0] ^ T
	#y_i = A ^ i * y_0
	#Y = [y_n-1, y_n-2,...y_0]
	#y = y_n
	#Y p = y
	#p = solve(Y, y)
	A = deepcopy(A)
	assert (A.shape[0] == A.shape[1])
	order = A.shape[0]
	if verbose == 1:
		print("###A###"), print(A)
		print("~~~Started computing characteristic polynom using krylov~~~")
	ys = []
	y = numpy.zeros((order, 1)).astype(float)
	y[0, 0] = 1.
	if verbose >= 2:
		print("###y_0###"), print(y), print("~~~Computed y_0~~~")
		if interactive: pause()
	ys.append(y)
	for i in range(1, order + 1):
		temp = numpy.dot(A, y)
		if verbose >= 2:
			print_group(A, y, temp), print("~~~Computing y_%d~~~" % i)
			if interactive: pause()
		y = temp
		ys.append(y)
	Y = numpy.concatenate(list(reversed(ys[:-1])), axis = 1)
	if verbose >= 1:
		print_group(Y, y), print("~~~Composed system for p~~~")
	if verbose >= 2:
		if interactive: pause()
	p = solver(Y, y, verbose = verbose, interactive = interactive)
	if verbose >= 1:
		print("###p###"), print(p)
		print("~~~Finshed computing characteristic polynom using krylov~~~")
	return numpy.poly1d(numpy.concatenate(([1.], -p.T[0])))
	return numpy.poly1d([1] + (-p.T[0]))

def polynom_danilevsky(A, *, verbose = 0, interactive = False):
	#To frobenius form
	A = deepcopy(A)
	assert(A.shape[0] == A.shape[1])
	order = A.shape[0]
	if verbose == 1:
		print("###A###"), print(A),
		print("~~~Started computing characteristic polynom using danilevsky~~~")
	S = numpy.identity(order)
	bound = order
	for k in range(1, order):
		B = numpy.identity(order)
		if abs(A[order - k, order -k - 1]) < 1e-9:
			pdb.set_trace()
			print("avoiding zero division")
			left_nonzero = numpy.nonzero(A[order-k, :order - k - 1])
			if len(left_nonzero[0]) == 0:
				#TODO: solve left top submatrix recursively
				pdb.set_trace()
				print("no left swapee")
				pass
			else:
				Q = numpy.identity(order)
				swap_cols(Q, left_nonzero[0][0], order - k - 1)
				if verbose >= 2:
					print_group(Q.T, B, Q)
					print("~~~Swapping rows %d with %d~~~" % 
							left_nonzero[0][0], order - k - 1
					)
					if interactive: pause()
		B[order - k - 1] = numpy.array([(1. if i == order-k-1 \
				else -A[order - k, i]) / A[order-k, order-k-1] for i in range(order)])
		S = numpy.dot(S, B)
		B_inv = numpy.identity(order)
		B_inv[order - k - 1] = A[order-k]
		if verbose >= 2:
			print_group(B_inv, A, B), print("~~~Doing step %d~~~" % k)
			if interactive: pause()
		A = numpy.dot(numpy.dot(B_inv, A), B)
	if verbose >= 2: print("###Composed B###"), print_group(S)
	if verbose >= 1:
		print("###A###"), print(A)
		print("~~~Finished computing characteristic polynom using danilevsky~~~")
	return numpy.poly1d(numpy.concatenate(([1.], -A[0])))

def derivative(f, eps=1e-6):
	return (lambda x: (f(x + eps / 2) - f(x - eps / 2)) / eps)

def root_newton(f, left, right, *, eps=1e-3, verbose = 0, interactive = False):
	if verbose >= 1:
		print("~~~Started newton method~~~")
		print("###left = %f, right = %f, eps = %f###" % (left, right, eps))
	x = (left + right) / 2.
	if (verbose >= 0): print("###x_0 = %f###" % x)
	k = 0
	target_eps = abs(eps)
	while True:
		x_next = x - f(x) / derivative(f, eps)(x)
		eps = abs(f(x_next))
		if verbose >= 2:
			print("x_%d = %e = %e - %e / %e (eps = %e)" %
					(k + 1, x_next, x, f(x), derivative(f, eps)(x), eps)
			)
			if interactive: pause()
		x = x_next
		k += 1
		if abs(eps) < target_eps:
			print("~~~Finished newton method~~~")
			return x

def build_spline(x, y, solver=solve_tridiagonal, verbose=0, interactive=False):
	"""Builds cubic spline by given points
	Algorithm:
		TODO
	"""
	x = numpy.array(x).astype(float).flatten()
	y = numpy.array(y).astype(float).flatten()
	assert(x.shape[0] == y.shape[0])
	n = x.shape[0]
	if verbose >= 1:
		fun_header()
		print_group(numpy.array([x, y]))
		if interactive: pause()
	h = numpy.array([x[i] - x[i-1] for i in range(1, n)])
	if verbose >= 1:
		print("~~~Done building h array~~~")
		print_group(numpy.array([h]))
		if interactive: pause()
	A = numpy.zeros((n, n)).astype(float)
	b = numpy.zeros((n, 1)).astype(float)
	A[0][0] = 1
	A[0][1] = 0
	b[0] = 0
	if verbose >= 2:
		print("A[0][0] = 1")
		print("A[0][1] = 0")
		print("b[0] = 0")
		if interactive: pause()
	A[n-1][n-1] = 1
	A[n-1][n-2] = 0
	b[n-1] = 0
	if verbose >= 2:
		print("A[%d][%d] = 1" % (n-1, n-1))
		print("A[%d][%d] = 0" % (n-1, n-2))
		print("b[%d] = 0" % (n-1))
		if interactive: pause()
	for i in range(1,n-1):
		A[i][i-1] = h[i-1] / 6.
		A[i][i] = (h[i-1] + h[i]) / 3.
		A[i][i+1] = h[i] / 6
		b[i] = (y[i+1] - y[i])/h[i] - (y[i] - y[i-1])/h[i-1]
		if verbose >= 2:
			print("A[%d][%d] = %f / 6. = %f" % (i, i-1, h[i-1], A[i][i-1]))
			print("A[%d][%d] = (%f + %f) / 3. = %f", (i, i, h[i-1], h[i], A[i][i]))
			print("A[%d][%d] = %f / 6. = %f" % (i, i+1, h[i], A[i][i+1]))
			print("b[%d] = (%f - %f) / %f - (%f - %f) / %f = %f" %
					(i, y[i+1], y[i], h[i], y[i], y[i-1], h[i-1], b[i])
			)
	if verbose >= 1:
		print("~~~Done building tridiagonal system")
		print("###(A|b)###")
		print_group(A, b)
		if interactive: pause()
	M = solver(A, b, verbose=verbose, interactive=interactive).flatten()
	if verbose >= 2:
		print("~~~System solved~~~")
		print("###M###")
		print_group(numpy.array([M]))
	k = numpy.zeros((n-1, 4)).astype(float)
	for i in range(n-1):
		k[i][3] = (M[i+1] - M[i]) / (6. * h[i])
		k[i][2] = M[i] / 2.
		k[i][0] = y[i]
		k[i][1] = (y[i+1] - y[i]) / h[i]  - (k[i][3] * h[i]**2 + k[i][2] * h[i])
		if verbose >= 2:
			print("k[%d][0] = %f" % (i, y[i]))
			print("k[%d][2] = %f" % (i, y[i]))
			print("k[%d][3] = (%f - %f) / (6. * %f) = %f" %
					(i, M[i+1], M[i], h[i], k[i][3])
			)
			print("k[%d][1] = (%f - %f) / %f - %f * %f ** 2 - %f * %f" %
					(i, y[i+1], y[i], h[i], k[i][3], h[i], k[i][2], h[i])
			)
			if interactive: pause()
	if verbose >= 1:
		print("~~~Done building spline~~~")
		print_group(k)
		if interactive: pause()
		print("~~~Finished building cubic spline~~~")
	return k

def get_interval(x, t):
	if t < x[0]: return 0
	for i in range(len(x)):
		if t < x[i]: return i-1
	return len(x) - 1

def get_cubic_value(k, x, t):
	return k[0] + k[1] * (t-x) + k[2] * (t-x)**2 + k[3] * (t-x)**3

def make_lambda_spline(x, k):
	x = x.flatten()[:-1]
	k = deepcopy(k)
	return lambda t: get_cubic_value(
			k[get_interval(x, t)], x[get_interval(x, t)], t
	)

def min_square_approx(x, y, n, solver=solve_gauss_noswap, verbose=0,
		interactive=False):
	"""Least squares interpolation"""
	x = x.flatten()
	y = y.flatten()
	assert(x.shape[0] == y.shape[0])
	N = x.shape[0]
	if verbose >= 1:
		print("fun_header")
	A = numpy.zeros((n, n))
	for k in range(n):
		for j in range(n):
			sum = 0.
			s = ""
			for i in range(N):
				sum += x[i]**k * x[i]**j
				s += "+ %f**%d * %f**%d " % (x[i], k, x[i], j)
			A[k][j] = sum
			if verbose >= 2:
				print("A[%d][%d] =%s %f" %
						(k, j, s + "= " if verbose >= 3 else "", A[k][j])
				)
				if interactive: pause()
	b = numpy.zeros((n, 1))
	for j in range(n):
		sum = 0.
		s = ""
		for i in range(N):
			sum += y[i] * x[i]**j
			s += "+ %f* %f**%d" % (y[i], x[i], j)
		b[j] = sum
		if verbose >= 2:
			print("b[%d] =%s %f" % (j, s + "= " if verbose >= 3 else "", b[j]))
	if verbose >= 1:
		print("~~~Done building System~~~")
		print("###A|b###")
		print_group(A,b)
		if interactive: pause()
	c = solver(A, b, verbose=verbose, interactive=interactive).flatten()
	if verbose >= 1:
		print("~~~System solved~~~")
		print("###c###")
		print_group(numpy.array([c]))
		print(numpy.poly1d(c))
		if interactive: pause()
		fun_footer()
	return c

def omega(x):
	result = numpy.poly1d([1.])
	for x_i in x:
		result *= numpy.poly1d([1, -x_i])
	return result

def div_diff(x, y):
	result = 0.
	for i in range(len(x)):
		part = float(y[i])
		for j in range(len(x)):
			if i != j:
				part /= x[i] - x[j]
		result += part
	return result

def div_diff_table_top(x, y, *, verbose=0, interactive=False):
	x = numpy.array(x).astype(float).flatten()
	y = numpy.array(y).astype(float).flatten()
	assert (x.shape == y.shape)
	n = len(x)
	f = numpy.zeros((n, n))
	f.fill(float("nan"))
	for i in range(n):
		f[i][0] = y[i]
		if verbose >= 1:
			print("###f[%d] = %f" % (i, f[i][0]))
			if interactive: pause()
	for j in range(1, n):
		for i in range(n-j):
			f[i][j] = (f[i+1][j-1] - f[i][j-1]) / (x[i+j] - x[i])
			if verbose >= 2:
				print("###f[%d--%d] = (%f - %f) / (%f - %f) = %f" %
						(i, i+j, f[i+1][j-1], f[i][j-1], x[i+j], x[i], f[i][j]) )
				if interactive: pause()
	return f
	if verbose >= 1:
		print("###f###")
		print_group(f)
		if interactive: pause()

def interpolate_newton(x, y, verbose = 0, interactive = False):
	"""Builds newton interpolation polynomial
	Algorithm:
		P_0(x) = y
		FUCK IT
		P_{n+1}(x) = P_{n}(x) + omega_{n}(x) * f_(x_0, ..., x_{n-1})
		omega_n(x) = Product_{k=0}^{n-1} (x - x_i)
	"""
	x = numpy.array(x).astype(float).flatten()
	y = numpy.array(y).astype(float).flatten()
	assert(x.shape[0] == y.shape[0])
	if verbose >= 1:
		fun_header()
		print("x---y")
		print_group(numpy.array([x, y]))
		f = div_diff_table_top(x, y, verbose=verbose, interactive=interactive)
	n = x.shape[0]
	result = numpy.poly1d([0.])
	for i in range(n):
		part = omega(x[:i]) * div_diff(x[:i+1], y[:i+1])
		if verbose >= 2:
			print("~~~Building item %d~~~" % i)
			print("###item###")
			print(part)
		result += part
	if verbose >= 1:
		print("~~~Done building polynomial~~~")
		print("###result###")
		print(result)
		fun_footer()
	return result

def fin_diff(y):
	if len(y) == 1:
		return y[0]
	else:
		return fin_diff(y[1:]) - fin_diff(y[0:-1])

def interpolate_newton_equal(x, y, verbose = 0, interactive = False):
	"""Builds newton interpolating polynomial assuming equal intervals"""
	x = numpy.array(x).astype(float).flatten()
	y = numpy.array(y).astype(float).flatten()
	#pdb.set_trace()
	assert (x.shape[0] == 2)
	n = y.shape[0]
	h = (x[-1] - x[0]) / (n-1)
	if verbose >= 1:
		fun_header()
		print("###x_start, x_end, n###")
		print(x[0], x[-1], n)
	result = numpy.poly1d([0.])
	omega = numpy.poly1d([1.])
	result += numpy.poly1d(y[0])
	for i in range(1, n-1):
		omega *= numpy.poly1d([1, -(i-1)])
		result += omega * fin_diff(y[:(i+1)]) / factorial(i)
	return result

def quadrature_trapezoid_points(x, y, verbose=0, interactive=False):
	"""Integrates using trapezoidal rule by given points"""
	x = numpy.array(x).astype(float).flatten()
	y = numpy.array(y).astype(float).flatten()
	assert(x.shape[0] == y.shape[0])
	if verbose >= 1:
		fun_header()
		print("###x###"), print(x)
		print("###y###"), print(y)
	n = x.shape[0]
	result = 0.0
	for i in range(n-1):
		if verbose >= 2:
			print("f_%d * dx_%d = %f * %f = %f" % (i, i, (y[i] + y[i+1]) / 2.,
					x[i+1] - x[i], (y[i] + y[i+1]) / 2. * (x[i+1] - x[i])
			))
			if interactive: pause()
		result += (y[i] + y[i+1]) / 2. * (x[i+1] - x[i])
	if verbose >= 1:
		print("###Integral = %f" % result)
		fun_footer()
	return result

def quadrature_function(a, b, f, n):
	x = [a + (b-a) * i / float(n+1) for i in range(n+1)]
	y = [f(t) for t in x]
	return quadrature_trapezoid_points(x, y)

def quadrature_max_accuracy_quotients(a, b, n, wp, *,
		solver=solve_gauss_noswap, verbose = 0, interactive = False):
	"""
	a, b --- left, right bounds
	n --- power of polynomial
	"""
	fun_header(), print("###a = %f, b= %f###, n = %d" % (a, b, n))
	a = float(a)
	b = float(b)
	A = numpy.zeros((n, n)).astype(float)
	g = numpy.zeros((n, 1)).astype(float)
	# wp --- additional weighting function
	for i in range(n): # ortogonalizing with x^i
		for j in range(n): # coefficient before x^j
			degree = i + j + wp + 1.
			A[i][j] = (b**degree - a**degree) / degree
		degree = i + n + wp + 1.
		g[i][0] = -(b**degree - a**degree) / degree
	if verbose >= 1:
		print("###A|g###"), print_group(A, g)
	x = solver(A, g, verbose=verbose, interactive=interactive)
	if verbose >= 1:
		print("###A|x|g###"), print_group(A, x, g)
	fun_footer()
	return x

def quadrature_max_accuracy(k, wp, *, verbose=0, interactive=True):
	"""
	"""
	fun_header()
	k = numpy.array(k).flatten().astype(float)
	n = k.shape[0]
	k = numpy.concatenate((k, [0]))[::-1]
	w = numpy.poly1d(k)
	wd = w.deriv()
	x = numpy.roots(w)
	A = numpy.zeros(n)
	weight = numpy.poly1d([1] + n * [0])
	for i in range(n):
		pass
	fun_footer()

def solve_integral_fredholm_singular(a, b, _lambda, alpha, beta, f, *,
		solver = solve_gauss_noswap, verbose = 0, interactive=True):
	assert (len(alpha) == len(beta))
	pdb.set_trace()
	if verbose >= 1: fun_header()
	n = len(alpha)
	B = numpy.zeros((n, n)).astype(float)
	g = numpy.zeros((n, 1)).astype(float)
	for i in range(n):
		for j in range(n):
			B[i, j] = quadrature_function(a, b,
					lambda s: beta[i](s) * alpha[j](s), 1000
			)
		g[i] = quadrature_function(a, b,
				lambda s: f(s) * beta[i](s), 1000
		)
	if verbose >= 1:
		print("###Beta|f###"), print_group(B, g)
		if interactive: pause()
	temp = numpy.identity(n) - _lambda * B
	A = solver(temp, _lambda * g, verbose=verbose, interactive=interactive)
	return A
	if verbose >= 1: fun_footer()

test.py

from solve_library import *
import numpy
import scipy
import matplotlib.pyplot as plot
from copy import deepcopy

def sum_squares(n):
	result = 0.
	for i in range(1, n+1):
		result += (2*i-1) ** 2
	return result

n = 7
s = numpy.zeros((n, n))
for i in range(n):
	s[i][0] = sum_squares(i+1)

for j in range(1, n):
	for i in range(j, n):
		s[i][j] = (s[i][j-1] - s[i-1][j-1]) / (j)

print(s)

units.py

from solve_library import *
import numpy
import scipy
import matplotlib.pyplot as plot
from copy import deepcopy

numpy.set_printoptions(precision = 4, suppress = True)
#numpy.set_printoptions(formatter={'float': '{: 0.2f}'.format})

print("Begin test print_group")
A = numpy.arange(1, 3*3+1).reshape(3, 3)
B = numpy.eye(3)
b = numpy.arange(2, 3+2).reshape(3, 1)
print("###A###"), print(A)
print("###B###"), print(B)
print("###b###"), print(b)
allinone = [A, B, b]
print("print_group(*allinone)"), print_group(*allinone)
print("print_group()"), print_group()
print("print_group(A, b)"), print_group(A, b)
print("print_group(A, B, b)"), print_group(A, B, b)
print("End test print_group")

print("Begin test swap_rows(M, i, j)")
M = numpy.arange(4 * 4).reshape(4, 4)
print("###M###"), print(M)
print("swap_rows(M, 2, 0)"), swap_rows(M, 2, 0)
print("###M###"), print(M)
print("End test swap_rows(M, i, j)")

print("Begin test swap_cols(M, i, j)")
M = numpy.arange(4 * 4).reshape(4, 4)
print("###M###"), print(M)
print("swap_cols(M, 2, 1)"), swap_cols(M, 2, 1)
print("###M###"), print(M)
print("End test swap_rows(M, i, j)")

print("Begin test to_triangle_top no_swap")
A = numpy.array([
	[2, 1, 1, 0],
	[4, 3, 3, 1],
	[8, 7, 9, 5],
	[6, 7, 9, 8],
])
print("###A###"), print_group(A)
P, Q = to_triangle_top(A, do_swap_rows = True, do_swap_cols = True, verbose = 2)
print("###P|A|Q|U###")
U = numpy.dot(numpy.dot(P, A), Q)
print_group(P, A, Q, U)
assert (numpy.linalg.norm(numpy.tril(U, -1) < 1e-6))
assert (numpy.linalg.norm(numpy.diag(U) > 1e-6))
print("End test to_triangle_top no_swap")

print("Begin test solve_triangle_top")
A = numpy.array([
	[3, 2, 1],
	[0, 3, 2],
	[0, 0, 3],
]).astype(float)
Q = numpy.array([
	[0, 0, 1],
	[0, 1, 0],
	[1, 0, 0],
]).astype(float)
b = numpy.dot(A, numpy.dot(Q.T, numpy.arange(1, 3 + 1).reshape(3, 1)))
print("###Q###"), print_group(Q)
x = solve_triangle_top(A, Q, b, verbose = 2)
print("###x###"), print_group(x)
assert (numpy.linalg.norm(x - numpy.arange(1, 3 + 1).reshape(3, 1)) < 1e-6)
print("End test solve_triangle_top")

print("Begin solve_triangle_bottom")
A = numpy.array([
	[3, 0, 0],
	[2, 3, 0],
	[1, 2, 3],
]).astype(float)
Q = numpy.array([
	[0, 0, 1],
	[1, 0, 0],
	[0, 1, 0],
]).astype(float)
b = numpy.dot(A, numpy.dot(Q.T, numpy.arange(1, 3 + 1).reshape(3, 1)))
print("###Q###"), print_group(Q)
x = solve_triangle_bottom(A, Q, b, verbose = 2)
print("###x###"), print_group(x)
assert (numpy.linalg.norm(x - numpy.arange(1, 3 + 1).reshape(3, 1)) < 1e-6)
print("End   solve_triangle_bottom")

print("Begin test solve_gauss")
A = numpy.array([
	[2, 1, 1, 0],
	[4, 3, 3, 1],
	[8, 7, 9, 5],
	[6, 7, 9, 8],
])
b = numpy.dot(A, numpy.arange(1, 4 + 1).reshape(4, 1))
print("###A|b###"), print_group(A, b)
x = solve_gauss(A, b, do_swap_rows=True, do_swap_cols=True, verbose = 2)
print("###x###"), print_group(x)
assert (numpy.linalg.norm(x - numpy.arange(1, 4 + 1).reshape(4, 1)) < 1e-6)
print("End test sovle_gauss")

#Archived
#print("Begin lu_decompose(A, verbose = 2)")
#A = numpy.array([
#	[2, 1, 1, 0],
#	[4, 3, 3, 1],
#	[8, 7, 9, 5],
#	[6, 7, 9, 8],
#]).astype(float)
#print("###A###"), print(A)
#L = lu_decompose(A, verbose = 2)
#print("###L###"), print(L)
#print("###A###"), print(A)
#print("End   lu_decompose(A, verbose = 2)")	

print("Begin lu_decompose(A, verbose = 2)")
A = numpy.dot(
	numpy.tril(numpy.arange(1, 4 * 4 + 1).reshape(4, 4), -1) + numpy.identity(4),
	numpy.triu(numpy.arange(1, 4 * 4 + 1).reshape(4, 4))
)
print("###A###"), print(A)
P, Q, L, U = lu_decompose(A, verbose = 2)
print("###P###"), print(P)
print("###A###"), print(A)
print("###Q###"), print(Q)
print("###L###"), print(L)
print("###U###"), print(U)
assert (numpy.linalg.norm(numpy.dot(L, U) - numpy.dot(numpy.dot(P, A), Q)) < 1e-6)
print("End   lu_decompose(A, verbose = 2)")

print("Begin solve_lu(A, b, verbose = 2)")
A = numpy.array([
	[2, 1, 1, 0],
	[4, 3, 3, 1],
	[8, 7, 9, 5],
	[6, 7, 9, 8],
]).astype(float)
b = numpy.dot(A, numpy.arange(1, 5).reshape(4, 1))
print("###b###"), print(b)
x = solve_lu(A, b, verbose = 2)
print("###x###"), print(x)
assert (numpy.linalg.norm(x - numpy.arange(1, 5).reshape(4, 1)) < 1e-6)
print("End    solve_lu(A, b, verbose = 2)")

###UPDATES####

print("Begin sqrt_decompose(A, verbose = 2)")
L = numpy.array([
	[1, 0, 0, 0],
	[1, 3, 0, 0],
	[1, 2, 3, 0],
	[1, 1, 1, 1],
]).astype(float)
A = numpy.dot(L, L.T)
print("###A###"), print(A)
L_computed = sqrt_decompose(A, verbose = 2)
print("###L|L^T|L*(L^T)|A###"), print_group(L_computed, L_computed, numpy.dot(L, L.T), A)
assert (numpy.linalg.norm(L_computed - L) < 1e-6)
print("End   sqrt_decompose(A, verbose = 2)")

print("Begin solve_sqrt(A, verbose = 2)")
L = numpy.array([
	[1, 0, 0, 0],
	[1, 3, 0, 0],
	[1, 2, 3, 0],
	[1, 1, 1, 1],
]).astype(float)
A = numpy.dot(L, L.T)
b = numpy.dot(A, numpy.arange(1, 4 + 1).reshape(4, 1))
x = solve_sqrt(A, b, verbose = 2)
print("###A * x = b###"), print_group(A, x, b)
assert (numpy.linalg.norm(x - numpy.arange(1, 4 + 1).reshape(4, 1)) < 1e-6)
print("End   solve_sqrt(A, verbose = 2)")

print("Begin solve_tridiagonal(A, b, verbose=2)")
A = numpy.array([
	[2, 1, 0, 0],
	[1, 10, -5, 0],
	[0, 1, -5, 2],
	[0, 0, 1, 4],
]).astype(float)
x_precise = numpy.array([-3, 1, 5, -8]).reshape(4, 1)
b = numpy.dot(A, x_precise)
print("A*x=b")
print_group(A, x_precise, b)
x = solve_tridiagonal(A, b, verbose=2, )
assert (numpy.linalg.norm(x - x_precise) < 1e-6)
print("End solve_tridiagonal(A, b, verbose=2)")

eps = 1e-3
weight = 0.75
print("Begin solve_jacobi(A, b, eps = %f, weight = %f, verbose = 2" % (eps, weight))
A = numpy.ones((4, 4)) + numpy.diag(numpy.arange(3, 4 + 3))
b = numpy.dot(A, numpy.arange(1, 4 + 1).reshape(4, 1))
print("###A|b###"), print_group(A, b)
x = solve_jacobi(A, b, eps = eps, weight = weight, verbose = 2)
print("###x###")
print(x)
assert (numpy.linalg.norm(x - numpy.arange(1, 4 + 1).reshape(4, 1)) < eps)
print("End   solve_jacobi(A, b, eps = %f verbose = 2)" % eps)

eps = 1e-3
print("Begin solve_seidel(A, b, eps = %f, verbose = 2" % eps)
A = numpy.ones((4, 4)) + numpy.diag(numpy.arange(3, 4 + 3))
b = numpy.dot(A, numpy.arange(1, 4 + 1).reshape(4, 1))
print("###A|b###"), print_group(A, b)
x = solve_seidel(A, b, eps = eps, verbose = 2)
print("###x###")
print(x)
assert (numpy.linalg.norm(x - numpy.arange(1, 4 + 1).reshape(4, 1)) < eps)
print("End   solve_seidel(A, b, eps = %f, verbose = 2)" % eps)

eps = 1e-3
print("Begin solve_gauss_seidel(A, b, eps = %f, verbose = 2" % eps)
A = numpy.ones((4, 4)) + numpy.diag(numpy.arange(3, 4 + 3))
b = numpy.dot(A, numpy.arange(1, 4 + 1).reshape(4, 1))
print("###A|b###"), print_group(A, b)
x = solve_gauss_seidel(A, b, eps = eps, verbose = 2)
print("###x###")
print(x)
assert (numpy.linalg.norm(x - numpy.arange(1, 4 + 1).reshape(4, 1)) < eps)
print("End   solve_gauss_seidel(A, b, eps = %f, verbose = 2" % eps)

eps = 1e-3
print("Begin solve_grad_decent(A, b, eps = %f, verbose = 2" % eps)
A = numpy.ones((4, 4)) + numpy.diag(numpy.arange(3, 4 + 3))
b = numpy.dot(A, numpy.arange(1, 4 + 1).reshape(4, 1))
print("###A|b###"), print_group(A, b)
x = solve_grad_decent(A, b, eps = eps, verbose = 2)
print("###x###")
print(x)
assert (numpy.linalg.norm(x - numpy.arange(1, 4 + 1).reshape(4, 1)) < eps)
print("End   solve_grad_decent(A, b, eps = %f, verbose = 2" % eps)

eps = 1e-3
print("Begin solve_coordinate_decent(A, b, eps = %f, verbose = 2" % eps)
A = numpy.ones((4, 4)) + numpy.diag(numpy.arange(3, 4 + 3))
b = numpy.dot(A, numpy.arange(1, 4 + 1).reshape(4, 1))
print("###A|b###"), print_group(A, b)
x = solve_coordinate_decent(A, b, eps = eps, verbose = 2)
print("###x###")
print(x)
assert (numpy.linalg.norm(x - numpy.arange(1, 4 + 1).reshape(4, 1)) < eps)
print("End   solve_coordinate_decent(A, b, eps = %f, verbose = 2" % eps)

############################ EXAM IS COMMING ############################

print("Begin rotation_matrix(order, i, j, theta)")
order, i, j, phi = 4, 0, 2, pi / 6
print("###order = %d, i = %d, j = %d, phi = %f###")
print("###rotation_matrix(order, i, j, theta)###")
print(rotation_matrix(order, i, j, phi))
print("End   rotation_matrix(order, i, j, theta)")

eps = 1e-3
print("Begin eig_rotation(A, eps = %f, verbose = 2)" % eps)
A = numpy.array([
	[ 10,  4, 10],
	[  4, -2,-14],
	[ 10,-14,  1]
])
print("###A###"), print(A)
print("###Spectrum###"),
v, H = eig_rotation(A, eps = eps, verbose = 3)
print("###Values###"), print(v)
print("###Vectors###"), print_group(H)
print("###eps=%f|obtained eps=%f###" % (eps, numpy.linalg.norm(numpy.dot(A, H) - v * H)))
assert (numpy.linalg.norm(numpy.dot(A, H) - v * H) < eps)
print("End   eig_rotation(A, eps = %f, verbose = 2)" % eps)

eps = 1e-3
print("Begin eig_max_power(A, eps = %f, verbose = 2)" % eps)
A = numpy.array([
	[ 10,  4, 10],
	[  4, -2,-14],
	[ 10,-14,  1]
])
print("###A###"), print(A)
eps = 1e-3
v, H = eig_max_power(A, eps = eps, verbose = 2)
print("###values###"), print(v)
print("###vectors###"), print(H)
print("End   eig_max_power(A, eps = %f, verbose = 2)" % eps)

print("Begin polynom_krylov(A, verbose = 2)")
A = numpy.array([
	[2.2, 1.0, 0.5, 2.0],
	[1.0, 1.3, 2.0, 1.0],
	[0.5, 2.0, 0.5, 1.6],
	[2.0, 1.0, 1.6, 2.0],
]).astype(float)
print("###A###"), print(A)
print(polynom_krylov(A, verbose = 2))
print("End   polynom_krylov(A, verbose = 2)")

print("Begin polynom_danilevsky(A, verbose = 2)")
A = numpy.array([
	[2.2, 1.0, 0.5, 2.0],
	[1.0, 1.3, 2.0, 1.0],
	[0.5, 2.0, 0.5, 1.6],
	[2.0, 1.0, 1.6, 2.0],
]).astype(float)
print("###A###"), print(A)
print(polynom_danilevsky(A, verbose = 2))
print("End   polynom_danilevsky(A, verbose = 2)")
#End per method unit tests

print("Begin derrivative(f, eps)")
f = lambda x: x * x - 3
for i in range(-1, 1 + 1):
	print("f'(%f) = %f" % (float(i), derivative(f)(float(i))))
print("End   derrivative(f, eps)")

print("Begin solve_newton(f, eps)")
f = lambda x: x * x - 3
result = root_newton(f, 1, 3, eps = 1e-3, verbose = 2)
print(result)
print("###error=%e###" % abs(result - sqrt(3)))
assert ( abs(result - sqrt(3)) < eps)
print("End solve_newton(f, eps)")

print("Begin test build_spline(x, y)")
x = numpy.array([0, 1, 3, 4, 6, 7])
y = numpy.array([0, 0, 1,-1, 0, 0])
k = build_spline(x, y, verbose=2)
print("###K###")
print_group(k)
fun = make_lambda_spline(x, k)
ix = numpy.linspace(0, 7, 8*10+1)
iy = numpy.array([fun(i) for i in ix])
plot.plot(x, y, 'o', ix, iy, '--')
#plot.show()
print("End build_spline(x, y)")

print("Begin min_square_approx(x, y, n)")
x = numpy.linspace(0.75, 3.75, 4+1)
y = numpy.array([2.50, 1.20, 1.12, 2.25, 4.28])
n = 2
c = min_square_approx(x, y, n, verbose=2)
p = numpy.poly1d(c[::-1])
print(c)
ix = numpy.linspace(0.75, 3.75, 4 * 10 + 1)
iy = p(ix)
plot.plot(x, y, 'o', ix, iy, '-')
#plot.show()
print("End   min_square_approx(x, y, n)")

print("Begin test div_diff_table_top(x, y)")
x = numpy.array([0, 1, 3]).astype(float)
y = numpy.array([2, 1, 4]).astype(float)
f = div_diff_table_top(x, y, verbose=2)
print("###f###")
print_group(f)
print("End test div_diff_table_top(x, y)")

print("Begin test interpolate_newton")
x = numpy.array([1, 2, 3, 5, 8])
y = numpy.array([1, -1, 1, -1, 1])
print("x---y")
print_group(numpy.array([x, y]))
p = interpolate_newton(x, y, verbose=2)
ix = numpy.linspace(x[0], x[-1], len(x)*10+1)
iy = p(ix)
plot.plot(x, y, 'o', ix, iy, '-')
#plot.show()
print("End test interpolate_newton")

print("Begin test interpolate_newton_equal")
x = [0, 2.5]
y = [1, 2, 1, 3, 1, 4]
h = (x[-1] - x[0]) / (len(y) - 1)
p = interpolate_newton_equal(x, y, verbose = 2)
x = numpy.linspace(0, 2.5, len(y))
f = lambda x: p((x - x[0]) / h)
ix = numpy.linspace(0, 2.5, (len(y) - 1) * 10 + 1)
iy = f(ix)
plot.plot(x, y, 'o', ix, iy, '-')
#plot.show()
print("End test interpolate_newton_equal")

print("Begin test quadrature_trapezoid_points")
x = [0, 0.3, 1.0]
y = [2, 1, 3]
integral = quadrature_trapezoid_points(x, y, verbose=2)
assert (abs(integral - 1.85) < 1e-6)
print("End   test quadrature_trapezoid_points")

print("Begin test quadrature_function")
print("End   test quadrature_function")

print("Begin quadrature_max_accuracy_quotients")
wp = 0
a = 0
b = 1
n = 3
k = quadrature_max_accuracy_quotients(a, b, n, wp, verbose=2)
result = numpy.array([-1/20., 3/5., -3/2.]).reshape((3, 1))
assert(numpy.linalg.norm(result - k) < 1e-6)
print("End   quadrature_max_accuracy_quotients")

print("Begin solve_integral_fredholm_singular")
a = 0.
b = 1.
alpha = [
	lambda x: x,
	lambda x: 1
]
beta = [
	lambda s: 1,
	lambda s: s
]
f = lambda x: 1 + x**3
_lambda = 1/2.
k = solve_integral_fredholm_singular(a, b, _lambda, alpha, beta, f, verbose=2, interactive=True)
print(k)
print("End   solve_integral_fredholm_singular")