import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import linprog
A_ub = np.array([[-2, -1],
[-2, -2],
[-1, -2],
[3, -1],
[-1, 3]])
b_ub = np.array([-2, -3, -2, 6, 6])
c = np.array([1, 1])
x1_bounds = (None, None)
x2_bounds = (None, None)
res = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=[x1_bounds, x2_bounds])
print(res)
con: array([], dtype=float64)
fun: 1.5000000000003295
message: 'Optimization terminated successfully.'
nit: 5
slack: array([2.50000000e-01, 6.59028387e-13, 2.50000000e-01, 4.50000000e+00,
4.50000000e+00])
status: 0
success: True
x: array([0.75, 0.75])
x = np.arange(1, 25) # 1, 2, ..., 24
y = np.array([75,77,76,73,69,68,63,59,57,55,54,52,50,50,49,49,49,50,54,56,59,63,67,72])
plt.figure(figsize=(7, 7))
plt.plot(x, y, 'o', label='Raw Data')
plt.xlabel('x', size=15)
plt.ylabel('y', size=15)
plt.title('Raw Data', size=18)
plt.grid()
plt.legend(loc='best', fontsize=15)
plt.show()
from scipy.optimize import fmin
from time import time
def f(a, x):
return a[0] * x**3 + a[1] * x**2 + a[2] * x + a[3]
def goal_function(a):
global x, y
return np.linalg.norm(f(a, x) - y)**2
a_start = np.array([1, 1, 1, 60])
ts = time()
results_fmin = fmin(goal_function, a_start)
time_fmin = time() - ts
Optimization terminated successfully.
Current function value: 37.492247
Iterations: 451
Function evaluations: 757
results_fmin
x = np.arange(1, 25) # 1, 2, ..., 24
y = np.array([75,77,76,73,69,68,63,59,57,55,54,52,50,50,49,49,49,50,54,56,59,63,67,72])
plt.figure(figsize=(7, 7))
plt.plot(x, y, 'o', label='Raw Data')
plt.plot(x, f(results_fmin, x), '-o', label='Fitted Data with fmin')
plt.xlabel('x', size=15)
plt.ylabel('y', size=15)
plt.title(r'Using $fmin$', size=18)
plt.grid()
plt.legend(loc='best', fontsize=15)
plt.show()
A = np.vstack([x**3, x**2, x, np.ones(x.shape)]).T
LSM_res = np.linalg.pinv(A) @ y
LSM_res
x = np.arange(1, 25) # 1, 2, ..., 24
y = np.array([75,77,76,73,69,68,63,59,57,55,54,52,50,50,49,49,49,50,54,56,59,63,67,72])
plt.figure(figsize=(7, 7))
plt.plot(x, y, 'o', label='Raw Data')
plt.plot(x, A @ LSM_res, '-', label='Fitted Data with LSM')
# plt.plot(x, f(results_fmin), 'o-', label='Fitted Data with fmin')
plt.xlabel('x', size=15)
plt.ylabel('y', size=15)
plt.title(r'Using LSM', size=18)
plt.grid()
plt.legend(loc='best', fontsize=15)
plt.show()
from scipy.optimize import differential_evolution
bounds = [(-100, 100), (-100, 100), (-100, 100), (-100, 100)]
a_start = np.array([1, 1, 1, 60])
x = np.arange(1, 25) # 1, 2, ..., 24
y = np.array([75,77,76,73,69,68,63,59,57,55,54,52,50,50,49,49,49,50,54,56,59,63,67,72])
def f(a, x):
return a[0] * x**3 + a[1] * x**2 + a[2] * x + a[3]
def goal_function(a):
global x, y
return np.linalg.norm(f(a, x) - y)**2
ts = time()
result_gen = differential_evolution(goal_function, bounds, x0=a_start,
strategy='currenttobest1bin',
mutation=[0.5, 1.8])
time_gen = time() - ts
print(result_gen)
fun: 37.492831874092936
jac: array([ 0.08455316, 0.00438405, -0.00508109, 0.04436916])
message: 'Optimization terminated successfully.'
nfev: 14125
nit: 233
success: True
x: array([ 9.17052893e-03, -1.58904090e-01, -1.75632012e+00, 8.02583608e+01])
plt.figure(figsize=(7, 7))
plt.plot(x, y, 'o', label='Raw Data')
plt.plot(x, f(result_gen.x, x), '-o', label='Fitted Data with genetic alg')
# plt.plot(x, f(results_fmin), 'o-', label='Fitted Data with fmin')
plt.xlabel('x', size=15)
plt.ylabel('y', size=15)
plt.title(r'Using Genetic Alg', size=18)
plt.grid()
plt.legend(loc='best', fontsize=15)
plt.show()
plt.figure(figsize=(7, 7))
plt.plot(x, y, 'o', label='Raw Data')
plt.plot(x, f(result_gen.x, x), '-', label='Fitted Data with genetic alg')
plt.plot(x, f(results_fmin, x), '-', label='Fitted Data with fmin')
plt.xlabel('x', size=15)
plt.ylabel('y', size=15)
plt.title(r'Comparing Genetic Alg and fmin', size=18)
plt.grid()
plt.legend(loc='best', fontsize=15)
plt.show()
plt.figure(figsize=(7, 7))
plt.plot(x, np.abs(f(result_gen.x, x) - f(results_fmin, x)), '-', label=r'$|f(a_{gen}) - f(a_{fmin})|$')
plt.xlabel('x', size=15)
plt.ylabel('y', size=15)
plt.title(r'Comparing Genetic Alg and fmin', size=18)
plt.grid()
plt.legend(loc='best', fontsize=15)
plt.show()
import pandas as pd
df = pd.DataFrame(data={'Genetic Alg':[time_gen, str(result_gen.nit)], 'fmin': [time_fmin, str(451)]}, index=['time', 'num iter'])
display(df)
from scipy.linalg import dft
np.set_printoptions(precision=2, suppress=True)
m = dft(8)
f_ = np.zeros((8, 1))
f_a = f_.copy()
f_a[:1] = 1
f_b = f_.copy()
f_b[:2] = 1
f_c = f_.copy()
f_c[:3] = 1
ff_a = m @ f_a
ff_b = m @ f_b
ff_c = m @ f_c
plt.figure(figsize=(7, 7))
plt.plot((ff_a * ff_a.conj()).real, '-o', label=r'$(a)$')
plt.plot((ff_b * ff_b.conj()).real, '-o', label=r'$(b)$')
plt.plot((ff_c * ff_c.conj()).real, '-o', label=r'$(c)$')
plt.title(r'$|\hat{f}|^2$', size=18)
plt.grid()
plt.legend(loc='best', fontsize=15)
plt.show()
def get_beta(k):
return 2 * (1 - (-1)**k) / (k * np.pi)
def get_b(k, omega2_0=2):
return get_beta(k) / (omega2_0 - k**2)
T = np.linspace(0, 10*np.pi, 500)
def x_a(t, omega2_0=2, N=9999):
res = 0.
for i in range(1, N):
k = i + 1
res += get_b(k) * (np.sin(k*t) - k / np.sqrt(omega2_0) * np.sin(np.sqrt(omega2_0) * t))
return res
def x_b(t, omega2_0=1, N=9999):
res = 0.
for i in range(1, N):
k = i + 1
if k == 1:
res += - get_beta(1) / 2 * t * np.cos(t)
res += get_beta(1) / 2 * np.sin(t)
else:
res += get_b(k, omega2_0=omega2_0) * (np.sin(k*t) - k / np.sqrt(omega2_0) * np.sin(np.sqrt(omega2_0) * t))
return res
x_a_res = [x_a(t) for t in T]
x_b_res = [x_b(t) for t in T]
plt.figure(figsize=(10, 5))
plt.plot(T, x_a_res, label=r'$omega_0^2=2$')
plt.plot(T, x_b_res, label=r'$omega_0^2=1$')
plt.title(r'$x(t)$', size=18)
plt.grid()
plt.legend(loc='best', fontsize=15)
plt.xlabel('t', size=15)
plt.show()
x = np.linspace(0, np.pi, 100)
def get_sin_thr_cos(x, N=10):
res = 2/np.pi
for i in range(N):
k = i + 1
if k!=1:
res += ((-1)**k + 1) * np.cos(k * x) / (np.pi * (1 - k**2)) * 2
return res
def get_cos_thr_sin(x, N=10):
res = 0.
for i in range(1, N):
k = i + 1
if k != 1:
res += k * ((-1)**k + 1) * np.sin(k * x) / (np.pi * (k**2 - 1)) * 2
return res
cos_ = get_cos_thr_sin(x)
sin_ = get_sin_thr_cos(x)
plt.figure(figsize=(10, 5))
plt.plot(x, sin_, label='Fourier series')
plt.plot(x, np.sin(x), label='true sin(x)')
plt.title('comparing sin(x)', size=18)
plt.grid()
plt.legend(loc='best', fontsize=15)
plt.xlabel('x', size=15)
plt.show()
plt.figure(figsize=(10, 5))
plt.plot(x, cos_, label='Fourier series')
plt.plot(x, np.cos(x), label='true cos(x)')
plt.title('comparing cos(x)', size=18)
plt.grid()
plt.legend(loc='best', fontsize=15)
plt.xlabel('x', size=15)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
!pip3 install --upgrade pip
!pip3 install --upgrade jax jaxlib
Requirement already satisfied: pip in /root/venv/lib/python3.7/site-packages (21.3.1)
Requirement already satisfied: jax in /root/venv/lib/python3.7/site-packages (0.2.25)
Requirement already satisfied: jaxlib in /root/venv/lib/python3.7/site-packages (0.1.74)
Requirement already satisfied: absl-py in /shared-libs/python3.7/py/lib/python3.7/site-packages (from jax) (0.15.0)
Requirement already satisfied: numpy>=1.18 in /shared-libs/python3.7/py/lib/python3.7/site-packages (from jax) (1.19.5)
Requirement already satisfied: typing-extensions in /shared-libs/python3.7/py-core/lib/python3.7/site-packages (from jax) (3.10.0.2)
Requirement already satisfied: opt-einsum in /shared-libs/python3.7/py/lib/python3.7/site-packages (from jax) (3.3.0)
Requirement already satisfied: scipy>=1.2.1 in /shared-libs/python3.7/py/lib/python3.7/site-packages (from jax) (1.7.2)
Requirement already satisfied: flatbuffers<3.0,>=1.12 in /shared-libs/python3.7/py/lib/python3.7/site-packages (from jaxlib) (1.12)
Requirement already satisfied: six in /shared-libs/python3.7/py-core/lib/python3.7/site-packages (from absl-py->jax) (1.16.0)
import jax.numpy as jnp
from tqdm.notebook import tqdm
xs = jnp.arange(1, 11)
ys = jnp.array([0.78,1.27,1.33,1.69,1.96,1.67,2.07,2.11,1.91,1.92])
def loss(x):
return jnp.power(ys - x[0] * (1 - jnp.exp(-x[1] * xs)), 2).sum()
def newton(x0, f, grad_f, hess_f, max_iter=100, eps=1e-6, verbose=False):
x_next = x0.copy()
ret_stat = {'fs': [],
'xs': [],
'grad_f_norm': []}
for i in range(max_iter):
h = jnp.linalg.solve(hess_f(x_next), grad_f(x_next))
x_next -= h
ret_stat['fs'].append(f(x_next))
ret_stat['xs'].append(x_next)
ret_stat['grad_f_norm'].append(jnp.linalg.norm(grad_f(x_next)))
if verbose:
print(f"a = {x_next[0]}, b = {x_next[1]}, f = ret_stat['fs'][-1]")
if ret_stat['grad_f_norm'][-1] < eps:
break
return ret_stat['xs'][-1], ret_stat['fs'][-1], ret_stat
import jax
from jax import grad, jacfwd, jacrev
def hessian(f):
return jacfwd(jacrev(f))
key = jax.random.PRNGKey(0)
#change x0 to obtain different solutions
x0 = jax.random.normal(key, (2,)) * 1#100, 5, 2, 3, 6, 10
res_x, res_f, ret_stat = newton(x0, loss, grad(loss), hessian(loss), max_iter=100, eps=1e-6, verbose=False)
res_x, res_f
x0 = jax.random.normal(key, (2,)) * 3 #100, 5, 2, 3, 6, 10
res_x, res_f, ret_stat = newton(x0, loss, grad(loss), hessian(loss), max_iter=100, eps=1e-6, verbose=False)
res_x, res_f
x0 = jax.random.normal(key, (2,)) * 5 #100, 5, 2, 3, 6, 10
res_x, res_f, ret_stat = newton(x0, loss, grad(loss), hessian(loss), max_iter=100, eps=1e-6, verbose=False)
res_x, res_f
x0 = jax.random.normal(key, (2,)) * 10 #100, 5, 2, 3, 6, 10
res_x, res_f, ret_stat = newton(x0, loss, grad(loss), hessian(loss), max_iter=100, eps=1e-6, verbose=False)
res_x, res_f