Using Sympy , Numpy, Scipy and Pandas for 3D - Plotting and Optimization¶
In [20]:
from sympy import init_printing
init_printing() #for pretty latex Output
import sympy as sp #Symbolic Math Library
import scipy as sc #Scientific Python Library for Multidimensional Optimization
from scipy.optimize import minimize #Importing Minimize
from sympy import * #Imports all special functions including cos, acos...
import numpy as np #Imports numeric python
import re #imports String manipulation and filtering functions
import pandas as pd #Importing for Dataframe - Excel Like Array - Printing
Plotting works only if there are 2 Independants (3D Plot)¶
In [3]:
x0, x1=sp.symbols('x0 x1') #define how many independent variables expression has!
In [4]:
F_string="Abs(((cos(x0))**4+(cos(x1))**4-2*((cos(x0))**2)*((cos(x1))**2))/(sqrt(x0**2+2*x1**2)))"
F_string
Out[4]:
In [5]:
F_sym=sp.sympify(F_string)
In [23]:
F_sym
Out[23]:
This Part will work if 2 independants were defined previously¶
In [7]:
F_lam=sp.lambdify((x0,x1),F_sym,modules='numpy')#Matplotlib needs a lambda function
x0_t=np.linspace(0.5,10)# sets ranges for a Plot
x1_t=np.linspace(0.5,10)
X0,X1=np.meshgrid(x0_t,x1_t)
Z=F_lam(X0,X1)
In [8]:
from mpl_toolkits.mplot3d.axes3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X0 , X1, Z, cmap='autumn', cstride=2, rstride=2)
plt.draw()
Building an Optimization Function from the String Function¶
Works with more than 2 Independants¶
In [9]:
#F_opt_str=F_string.replace("x0","x[0]").replace("x1","x[1]")
#Finds indices and converts to a scipy.optimize expression
indices = sorted(set(int(n) for n in re.findall(r'\bx(\d+)', F_string)),reverse=True)#reverse is important because otherwise x[10] becomes x[1]0
#Set takes unique objects (erases dublicates)
t = F_string[:]
#global space
space=len(indices)
for ind in indices:
t = t.replace('x' + str(ind), 'x[' + str(ind) + ']')# if variable begins with "1" lower second to "ind-1"
F_opt_str=t
In [10]:
F_opt_str
Out[10]:
In [11]:
F_opt_lam=eval('lambda x:'+F_opt_str)#establish a lambda function of string_function
In [12]:
space
Out[12]:
In [14]:
#allows to See Iteration Progress: Method "COBYLA" does not allow callback functions
class iter_progress:
x_it=[]#allowing to access objects "outside functions"
y_it=[]
iter_log=[]
def __init__(self):
pass#iter_progress.x_iter
def build_iter(self,space):
import numpy as np
iter_progress.y_it=np.zeros((1,1), dtype=float)
iter_progress.x_it=np.zeros((1,space), dtype=float)
iter_progress.iter_log=np.zeros((1,space+1), dtype=float)
def obj(self,x):
#x=np.array([1,2])
out=F_opt_lam(x)
# print(out)
out=np.array([[out]])
x_copy=x.copy()[None]
#print("x_it=")
iter_progress.x_it=np.concatenate((iter_progress.x_it, x_copy))
#print(iter_progress.x_it)
#print(" ")
iter_progress.y_it=np.concatenate((iter_progress.y_it, out))
#print(iter_progress.y_it)
#print(" ")
#print(out)
return out
def mine(self,x0=None):
if x0==None:# if no starting value is defined "a zero array is defined" does not work always
x0=np.zeros(((1,space)))
print("x0=%s" % x0)
#ip=iter_progress()
#ip.build_iter
res=minimize(ip.obj,x0=x0,method='SLSQP')
iter_progress.iter_log=np.concatenate((iter_progress.x_it, iter_progress.y_it),axis=1)
#print(iter_progress.x_it)
print(" ")
#print(iter_progress.y_it)
print(iter_progress.iter_log)
print(" ")
print(res.x)
In [15]:
ip=iter_progress()
ip.build_iter(space)
xt = np.array([[4,6]])
ip.mine(xt)#Define an optional starting x0 if you want - x0 = np.array([[1,2,...length space !]])
Simple Optimization without Using the Class¶
In [16]:
#Simple Optimization without Using the Class and Progress
#x0=np.array([3,4])
#res=minimize(F_opt_lam,x0=x0,method='SLSQP')
#res.x
#res
"Pandas" DataFrame for a "pretty" table-view¶
In [17]:
var_str=[]# strings of Columns for pandas
indices_help=indices[::1]#reverse indices
for ind in indices_help:
var_str.append('x'+str(ind))
var_str.append('f(x_i,n)')
In [18]:
var_str
Out[18]:
In [21]:
df=pd.DataFrame(iter_progress.iter_log, columns=var_str) #Pretty Version of Optimization
df
Out[21]:
In [22]:
from mpl_toolkits.mplot3d.axes3d import Axes3D
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(X0 , X1, Z, cmap='autumn', cstride=2, rstride=2)
ax.plot(iter_progress.x_it[1:,0],iter_progress.x_it[1:,1], zs=iter_progress.y_it[1:,0], zdir='z', label='path')
ax.view_init(20,50)
plt.axis([3.1,3.9,4.5,5.8])
plt.draw()
Keine Kommentare:
Kommentar veröffentlichen